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1. 


INTRODUCTION 


1 . 1  BACKGROUND 


The  detailed  structure  of  the  gravity  field  is  infer¬ 
red  from  data  collected  by  a  variety  of  instruments  (sensors). 
These  sensors  measure  different  quantities  and  have  different 
spectral  responses  and  noise  characteristics.  Gravity  esti¬ 
mates  are  obtained  by  combining  measurements  from  all  sources. 
From  a  practical  point  of  view,  it  is  of  fundamental  importance 


•  To  determine  the  accuracy  of  gravity 
estimates  obtained  from  data  already 
collected 

•  To  determine  what  additional  data  and  sur¬ 
vey  characteristics  would  yield  a  required 
accuracy  of  the  estimates  of  gravity  in 

a  given  region. 


1.2  STUDY  OBJECTIVES 

The  objectives  of  this  study  were  to  develop  a  method¬ 
ology  and  a  computer  program  for  quantifying  the  errors  in  the 
estimates  of  gravity  available  from  multisensor  survey  data. 
Survey  types  explicitly  considered  in  this  study  consist  of 
any  combination  of 

•  Satellite  radar  altimetry 

•  Satellite-to-satellite  tracking  (SST)  in 
a  high-low  configuration 
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•  Land-based/shipborne  gravimetry 

•  Airborne  gradiometry. 

Other  survey  possibilities  can  be  easily  incorporated  into  the 
formulation . 

A  separate  phase  of  this  study  dealt  with  the  estima¬ 
tion  of  weapons  systems'  target  miss  induced  by  errors  in  the 
gravity  estimates  obtained  from  multisensor  survey  data.  The 
results  of  this  second  phase  of  the  study  are  presented  in  a 
separate,  expanded  version  of  this  report. 


1.3  TECHNICAL  APPROACH 

The  approach  utilizes  a  flat-earth  approximation  that 
furnishes  algebraic  frequency-domain  interrelations  among  geo¬ 
detic  quantities.  An  extension  of  classical  Wiener  smoothing 
techniques  (Ref.  1)  permits  the  computation  of  the  average 
power  spectral  density  of  the  post-survey  gravity  residuals. 


Since  the  Wiener  smoother  is  optimal  in  the  sense  of 
minimizing  the  mean-square  residuals,  the  results  can  be  viewed 
as  representing  the  best  possible  use  of  the  survey  data. 
Actually,  there  are  some  approximations  involved  in  the  evalu¬ 
ation  of  the  spectral  density  of  the  residuals  and,  more  fun¬ 
damentally,  in  the  use  of  the  flat-earth  formulation  with  its 
entailing  loss  of  accuracy  for  wavelengths  comparable  to  the 
radius  of  the  earth.  Nonetheless,  since  most  of  the  energy  in 
the  gravity  anomaly,  deflections  of  the  vertical,  gravity  dis¬ 
turbance  vector  and  the  gradients  of  gravity  is  contained  in 
the  frequency  band  where  the  method  yields  accurate  results, 
the  analysis  does  provide  an  accurate  measure  of  the  perform- 
able  accuracy  of  the  estimates  of  these  quantities  obtained 
from  multisensor  survey  data. 


Figure  1.3-1  presents  a  graphic  illustration  of  the 
multisensor  survey  analysis  methodology.  The  circle  repre¬ 
sents  the  results  that  are  obtained  through  the  application  of 
the  techniques  discussed  herein.  They  are  the  statistics  of 
the  post-survey  residual  gravity  errors.  The  ovals  represent 
the  quantities  that  must  be  specified  in  the  evaluation  of  the 
errors.  These  are:  a  statistical  model  for  the  unsurveyed 
anomalous  field,  sensor  error  models,  and  the  characteristics 
of  the  survey.  The  box  represents  the  operations  performed  on 
the  input  specifications  to  obtain  the  statistics  of  the  resid¬ 
ual  gravity  field. 


R  62394a 


Figure  1.3-1  Multisensor  Survey  Analysis  Methodology 

The  field  model  is  the  a  priori  power  spectral  density 
of  the  anomalous  potential  at  the  earth's  surface.  The  statis¬ 
tics  of  all  other  field-related  quantities  are  obtained  from 
this  power  spectral  density  through  the  use  of  flat-earth  re¬ 
lations  (Ref.  2).  The  field  model  describes  the  local  behavior 
of  the  gravity  field  in  the  vicinity  of  the  region  of  interest 
(e.g.,  a  missile  launch  point).  Either  analytical  or  empirical 
models  can  be  used  in  the  analysis. 


I 


A  computer  program  was  designed  on  the  basis  of  the 
theory  discussed  in  this  report.  A  wide  variety  of  survey 
possibilities  were  simulated.  The  results  of  these  simula¬ 
tions  are  included  and  discussed  in  this  report. 


1.4  REPORT  ORGANIZATION 

The  organization  of  this  document  is  as  follows: 
Chapter  2  presents  an  overview  of  the  methodology  for  the  analy¬ 
sis  of  multisensor  gravity  survey  residual  errors  including 
error  models  for  the  various  survey  possibilities  considered. 
Chapter  3  presents  a  collection  of  results  obtained  through 
the  use  of  the  techniques  described  in  the  previous  chapter. 
Chapter  4  presents  a  summary  of  the  information  contained  in 
this  report  and  discusses  several  possible  extensions  of  the 
analysis . 


Various  appendices  complement  this  report.  Appendix  A 
presents  a  succinct  compilation  of  those  concepts  of  Fourier 
analysis  necessary  for  the  development  of  the  theory  of  Chap¬ 
ter  2.  Appendices  B  and  C  discuss  additional  technical  aspects 
of  the  analysis.  Appendix  B  gives  the  derivations  of  the  flat- 
earth  frequency-domain  relations.  Appendix  C  presents  detailed 
derivations  of  the  formulas  for  the  spectral  density  of  the 
post-survey  residuals. 


1.5  A  NOTE  ON  TERMINOLOGY 

It  is  customary  to  refer  to  functions  that  assume  ran¬ 
dom  values  as  random  processes  when  these  functions  depend  on  a 
single  variable  or  coordinate  and  as  random  fields  when  they 
depend  on  two  or  more  variables  or  coordinates  (see,  for  exam¬ 
ple,  Ref.  3).  To  avoid  confusion,  the  noun  field  is  used  only 
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in  reference  to  the  gravity  field.  Functions  that  assume  ran¬ 
dom  values  are  referred  to  as  random  processes  in  n-dimensional 
space  or  simply  as  random  processes  when  the  number  of  inde¬ 
pendent  variables  is  clear  from  the  context. 

Random  processes  can  be  scalar  or  vector-valued  depend 
ing  on  whether  to  every  point  in  space  there  corresponds  one 
or  more  than  one  real  numbers.  For  the  purpose  of  this  report, 
a  vector  denotes  any  list  of  numbers  associated  with  every 
point  of  the  coordinate  system.  This  list  of  numbers  need  not 
be  a  physical  vector.  For  example,  a  three-dimensional  vector 
process  in  two-dimensional  space  can  consist  of  the  undulation 
of  the  geoid,  the  gravity  anomaly  and  the  north  component  of 
the  gravity  disturbance  vector  over  the  earth's  surface. 


The  term  geodetic  quantities  refers  to  the  collection 
consisting  of  the  anomalous  potential,  the  undulation  of  the 
geoid,  the  deflections  of  the  vertical,  the  gravity  anomaly, 
the  gravity  disturbance,  the  components  of  the  gravity  dis¬ 
turbance  vector  and  the  gravity  gradients.  Linear  combina¬ 
tions  of  these  quantities  are  referred  to  as  field-related 
quantities . 


ANALYSIS  OF  MULTISENSOR  GRAVIMETRIC  SURVEY 
RESIDUAL  ERRORS 


This  chapter  presents  a  method  for  determining  the 
statistics  of  the  errors  in  the  anomalous  field  estimates  ob¬ 
tained  from  multisensor  survey  data.  The  method,  which  is 
based  on  an  extension  of  optimal  Wiener  smoothing  theory,  per¬ 
mits  a  characterization  of  the  errors  in  estimates  of  point 
and  spatial  averages  of  the  gravity  field  in  wavelengths  which 
are  small  compared  to  the  radius  of  the  earth.  Survey  types 
include  any  combination  of  satellite  altimetry,  satell ite- to- 
satellite  tracking  (SST) ,  land-based/shipborne  gravimetry  and 
airborne  gradiometry. 

In  the  analysis,  the  round  earth  is  approximated  by 
an  infinite  plane  (flat-earth)  which  practically  coincides 
with  the  earth  in  the  neighborhood  of  the  region  in  which  esti¬ 
mates  are  sought.  Frequency-domain  mappings  (transfer  func¬ 
tions)  are  used  to  related  field-related  variables  in  a  concise 
algebraic  manner.  The  anomalous  gravity  field  is  viewed  as  a 
realization  of  a  stationary  random  process  on  the  earth  plane. 
Values  of  the  field  above  the  earth  are  obtained  via  the  flat- 
earth  upward  continuation  formula  of  Heiskanen  and  Moritz 
(Ref.  2). 


Frequency-domain  techniques  are  extensively  used  in 
this  chapter.  Appendix  A  contains  a  summary  of  the  relevant 
concepts,  definitions,  and  results  of  Fourier  analysis. 

This  chapter  is  organized  as  follows:  Section  2 . 1 
presents  the  frequency-domain  relations  between  geodetic  quan¬ 
tities.  The  expression  for  the  average  PSD  of  the  post-survey 
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residuals  is  discussed  in  Section  2.2.  Survey  geometry  and 
error  models  are  described  in  Section  2.3. 

2.1  FREQUENCY-DOMAIN  FLAT-EARTH  RELATIONS 

To  determine  the  residual  gravity  errors  resulting 
from  the  performance  of  a  multisensor  survey,  it  is  necessary 
to  relate  in  a  geophysically  consistent  manner  the  different 
types  of  gravity  data  collected  by  the  survey. 

Geodetic  quantities  are  interrelated  by  integro- 
differential  operators.  For  example,  the  deflections  of  the 
vertical  are  given  in  terms  of  the  gravity  anomaly  by  the  for¬ 
mulas  of  Vening  Meinesz,  the  gravity  disturbance  vector  is  the 
gradient  of  the  anomalous  potential,  etc.  These  relations 
imply  specific  correlation  structures  when  the  gravity  field 
is  viewed  as  a  realization  of  a  stationary  random  field. 

A  similar  observation  can  be  made  for  the  post-survey 
residual  errors.  Since  the  errors  in  the  estimates  obtained 
with  any  geophysically  consistent  data  processing  algorithm  must 
satisfy  the  same  mathematical  constraints  linking  the  estimated 
variables,  the  statistics  of  the  errors  will  display  the  same 
correlation  structure  present  in  the  original  quantities. 

With  the  use  of  the  flat-earth  approximation,  the 
relatively  complex  integro-di f ferential  operators  become  simple 
algebraic  relations  between  Fourier  transforms.  Ratios  between 
Fourier  transforms  of  geodetic  quantities  turn  out  to  be  ra- 

JL. 

tional  functions  of  frequency  denominated  transfer  functions'. 

In  the  case  in  which  the  field  is  seen  as  a  realization  of  a 

*The  term  "transfer  function"  is  borrowed  from  an  analogous 
concept  in  Systems  Theory. 


random  process,  ratios  and  products  of  transfer  functions  re¬ 
late  the  power  spectral  and  cross-spectral  densities  of  all 
field-related  quantities. 

Table  2.1-1  presents  in  the  last  column  the  transfer 
functions  from  anomalous  surface  potential,  T  ,  to  the  geo¬ 
detic  quantities  listed  in  the  first  column  of  the  table.  The 
notation  used  throughout  this  report  for  all  geodetic  quanti¬ 
ties  is  given  in  the  second  column  of  the  table,  while  the 
third  column  presents  the  space-domain  relations  between  each 
of  the  variables  and  the  anomalous  potential.  The  constant  gQ 

used  in  the  table  represents  the  mean  value  of  gravity  of  the 

o 

earth  (gQ  -  9.798  m/sec  ).  The  interpretation  of  each  of  the 
transfer  functions  in  Table  2.1-1  is  discussed  next. 

/s 

Let  Tq(s)  be  the  Fourier  transform  of  the  anomalous 
surface  potential,  i.e., 

00 

Tq(s)  =  To(x)e”i27I<-’— >  dxxdx2  (2.1-1) 

-00 

T 

In  the  above  equation,  x  =  (x^,x2)  denotes  the  position  of  a 
point  on  the  plane  with  x^  and  x2  measured  in  the  east  and 
north  directions,  respectively,  s  =  (s^,s2)^  is  the  vector 
planar  frequency  with  s^  and  s2  representing  frequencies  meas¬ 
ured  in  the  east  and  north  directions,  i  =  /-T,  and  <x,s>  is 
the  inner  product  of  the  vectors  x  and  s. 

The  Fourier  transform  of  any  of  the  quantities  in  the 
first  column  of  Table  2.1-1  can  be  obtained  by  multiplying 
Tq(s)  by  the  corresponding  transfer  function  in  the  last  column 
of  the  table.  For  example,  the  Fourier  transform  of  the  anoma¬ 
lous  potential  at  height  z  is  given  by 


2-3 


TABLE  2.1-1 

TRANSFER  FUNCTIONS  FROM  ANOMALOUS  SURFACE  POTENTIAL 


. . — . - 

QUANTITY 

SYMBOL 

RELATION  TO  ANOMALOUS 
POTENTIAL 

TRANSFER  FUNCTION 
FROM  ANOMALOUS  * 
SURFACE  POTENTIAL 

Anomalous  Potential  at  Height  z 

Tz 

Tz 

e*2nsz 

Undulation  of  the  Geoid 

N 

V*o 

l/*o 

East  Deflection  of  the  Vertical 

n 

■(3T0/ax1)/g0 

-i2nsi/g0 

North  Deflection  of  the  Vertical 

i 

-OT0/3x2)/g0 

-i2ns2/go 

Gravity  Anomaly 

Ag 

-OTz/az)lz=0  -  2Tq/R 

2ns  -  2/R 

East  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

H 

aV3xl 

.  -2nsz 

i2ns^e 

North  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

1 

3Tz/3x2 

i2ns2e'2nsz 

Vertical  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

HH 

3Tz/3z 

-2nse’2nsz 

Gravity  Disturbance  at  Height  z 

-9T  /3z 

2nse"2nsz 

East-East  Gradient  at  Height  z 

2  2 

3  Tz/3xl 

,  2  2  -2nsz 

Sje 

North-North  Gradient  at  Height  z 

r*2x2 

2  2 

3*Tz/3x2 

-4n2s2e~2nsz 

Vertical-Vertical  Gradient  at 

Height  z 

rzz 

2  2 

3^Tz/3z* 

4n2s2e~2,ISZ 

East-North  Gradient  at  Height  z 

Pxlx2 

2 

3  Tz/3x13X2 

,  2  -2nsz 

s^S2e 

East-Vertical  Gradient  at  Height  z 

rXjZ 

32Tz/3x13z 

./  2  „  -2nsz 

-imh  s^se 

North-Vertical  Gradient  at  Height  z 

r  , 

A  2  “ 

32Tz/8x232 

./  2  -2nsz 

S2se 

T  ( s )  =  e'2nsz  T  ( s )  (2.1-2) 

A  complete  derivation  of  the  transfer  functions  presented  in 
Table  2.1-1  is  given  in  Appendix  B.l.  In  particular,  it  is 
shown  there  that  Eq.  2.1-2  is  the  frequency-domain  equivalent 
of  the  flat-earth  upward  continuation  formula  of  Ref.  2. 


2-4 


In  the  case  in  which  the  surface  anomalous  potential 
is  seen  as  a  realization  of  a  stationary  random  process,  all 
field-related  quantities  turn  out  to  be  jointly  stationary. 
In  fact,  it  is  shown  in  Appendix  B.2  that  if  w  and  u  represent 
two  vectors  of  field-related  quantities  with 

w  =  (w^  w2 . wp)T  (2.1-3) 

and 

u  =  (u-p  u2,  ....  uq)T  (2.1-4) 


then  the  cross-spectral  density  matrix  of  w  and  u,  4>  (s),  is 

w ,  u 

given  in  terms  of  the  power  spectral  density  (PSD)  oT  the  sur¬ 
face  anomalous  potential,  4>T  T  (s) ,  by^ 

o’  lo 


<P 

w,u 


_  jc 

G  F  <J>T  T 
1  o  * 1  o 


(2.1-5) 


where 

G(s)  =  [G1(s),  G2(s),  Gp(s) JT  (2.1-6) 

and 

F(s)  =  [F1(s),  F2(s),  Fq(s)]T  (2.1-7) 

are  the  vector  transfer  functions  from  surface  anomalous  poten¬ 
tial  to  w  and  u,  respectively. 

As  an  example,  consider  the  situation  where  the  vector 
w  consists  of  three  quantities:  the  east  deflection  of  the 
vertical,  the  north  deflection  of  the  vertical  and  the  undula¬ 
tion  of  the  geoid  so  that 

fA  superscript  asterisk  denotes  complex  conjugate  when  attached 
to  a  scalar  quantity  and  conjugate  transpose  when  attached  to 
a  matrix  or  a  vector. 
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(2.1-8) 


w  =  (n ,  i,  n)l 

and  suppose  the  vector  u  has  two  components:  the  undulation 
of  the  geoid  and  the  gravity  disturbance  at  height  z;  i.e., 

u  =  (N,  6g)T  (2.1-9) 

The  vector  transfer  functions  G  and  F  are,  from  Table  2.1-1, 


G(s)  =  (-i2ns1/gQ,  -i27ts2/go,  l/gQ) 


(2.1-10) 


and 


F(s)  =  (l/go,  2nse'27lsz)T 


Therefore,  from  Eq.  2.1-5,  the  cross-spectral  density  matrix 


<t> 


w,u 


^Vn 

*e,sg 

\  *N,N 

***N  ,6g 

(2.1-11) 


m 


v 


is  given  in  terms  of  the  PSD  of  the  anomalous  surface  potential, 

t  i  by 
o  ’  Ao 


4>  ( s ) 

w ,u 7 


/-i2ns1/g 


o 


-i2ns0/g  4 
2'  6o 

1/g  2 


-i47T2s1se_27TSZ/go\ 

-i4n2s9se”2nsz/g 


~  „  -2nsz, 
2nse  /g 


T  (  ) 

Lo’  lo 


°  / 


(2.1-12) 
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2.2 


MULTISENSOR  SURVEY  ERROR  ANALYSIS 


The  frequency-domain  flat-earth  relations  discussed 
in  the  previous  section  together  with  an  extension  of  Wiener 
smoothing  techniques  form  the  basis  for  the  analysis  of  the 
performable  accuracy  of  the  gravity  estimates  obtained  from 
multisensor  survey  data. 

The  problem  of  determining  the  errors  in  the  estima¬ 
tion  of  the  gravity  field  from  multisensor  survey  data  can  be 
formulated  in  the  following  manner:  it  is  desired  to  charac¬ 
terize  the  differences  (estimation  errors) 

6w(x)  =  w(x)  -  w°(x)  (2.2-1) 

between  the  true  values  of  the  process  w(x) ,  and  the  best  esti- 
mates  w  (x) .  The  components  of  w  are  any  collection  of  field- 
related  quantities. 

The  estimates  w°(x)  are  optimally  obtained  from  n 
classes  of  measurements 


rp  =  f£p(x) |x  e  Mp} ;  p  =  1 ,2 ....  ,n 


(2.2-2) 


corresponding  to  the  measurements  that  constitute  the  survey. 

Each  class  ro  represents  a  collection  of  vector  measurements, 
P 

»^p ,  in  which  all  measurement  points  form  a  rectangular  grid 

Mp  . 


The  data,  ^ ,  are  linear  combinations  of  field-related 

quantities,  uD ,  corrupted  by  additive  noise,  E~  ;  i.e., 

-p  ~p 


*In  the  mean-square  sense. 


(2.2-3) 


!kp(x)  =  up(x)  +  Ep(x);  x  £  Mp 

In  addition,  the  measurement  errors,  EQ ,  0  =  l,2,...,n,  are 

~ I p 

taken  to  be  stationary  gaussian  processes  independent  of  the 
gravity  field. 


The  grid  (see  Fig.  2.2-1)  is  completely  character¬ 
ized  by  three  quantities:  a  rotation  matrix,  0D ,  determined 

P 

by  the  orientation  of  the  survey  tracks  with  respect  to  the 
east  direction,  a  translation  vector  r^  which  locates  the  posi¬ 
tion  of  the  origin  of  the  grid  in  the  east-north  frame,  and  a 
spacing  matrix  ,  determined  by  the  separation  between  survey 
tracks  and  between  samples  along  a  track.  In  terms  of  the 

quantities  in  Fig.  2.2-1,  0O ,  r_  and  JQ  are  defined  by 

P  “ p  P 


and 


(cos  0  -sin  0 

sin  0  cos  0 


(2.2-4) 

(2.2-5) 


(2.2-6) 


For 

data  can  be 


the  survey  possibilities  considered  in  this 

divided  into  several  classes,  ro.  These  are 

P 


report , 


Two  scalar-measurement  classes  for  each 
satellite  radar-altimeter  mission.  One 
class  for  the  ascending  passes  and  another 
class  for  the  descending  passes 


R -62449 


X2 

(NORTH) 


Figure  2.2-1  Definition  of  the  Grid  M 


•  Two  scalar-measurement  classes  for  a 
high-low  satellite-to-satellite  tracking 
mission.  One  class  for  the  ascending 
and  another  for  the  descending  passes  of 
the  low  satellite 

•  One  scalar-measurement  class  for  a  land- 
based  or  ship  gravimetric  survey.  The 
differences  between  land-based  and  ship 
surveys  are  reflected  in  the  error  models 

•  One  vector-measurement  class  for  an  air¬ 
borne  gradiometer  survey.  Each  vector 
measurement  in  this  class  consists  of 
six  entries  corresponding  to  the  two 
outputs  from  each  of  the  instruments  in 
a  gradiometer  triad. 


Appendix  C  presents  a  detailed  analysis  of  the  problem 
formulated  above.  The  statistics  of  the  post-survey  residual 
errors  at  any  given  point  are  shown  to  depend  upon  the  relative 
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position  of  the  point  with  respect  to  the  measurements.  The 
average  spectral  density  of  the  residuals  is  then  defined  by 
averaging  over  all  possible  relative  locations  with  respect  to 
the  measurements.  For  the  residual  surface  anomalous  potential, 
the  average  spectral  density  is  given  by 


(f) 

V6T  ,6T 
o  *  o 


1 


+  1/4>t  t 
xo  ’  Lo 


(2.2-7) 


where  F„  is  the  vector  transfer  function  from  surface  anomalous 
~P 

potential  to  the  quantities  measured  in  class  p  ,  4><.  <.  is  the 

spectral  density  matrix  of  the  errors  in  class  p  (see  below) 

and  <J>T  T  is  the  PSD  of  the  unsurveyed  surface  anomalous  po- 
o  ’  o 

tential  obtained  from  a  gravity- field  model. 


The  spectral  density  matrix  of  the  errors  in  class  p 
consists  of  two  terms: 


<t> 


r  (s) 


=  ‘frp  F  (s) 

V-e  - 


*a  a 

£P’~P 


(2.2-8) 


The  first  term,  ,  is  the  unnormalized" 

-P’-P 

of  the  (discrete)  measurement  error  process 
Explicit  expressions  for  „  are  given  in 

-P’-P 

the  various  survey  possibilities  considered. 

<*>  _  *  is  the  spectral  density  matrix  of  the 

-p  ’-p 

in  class  p  computed  from 


spectral  density 

Eft  of  Eq.  2.2-3. 
Section  2.3  for 

The  second  term 
aliasing  errors 


1 


*  .  (s)  =  FR(s+eRj'1A)  Fr  ( s+6ft  JrXA  )  4>T  T  (s+eRj:1A) 

-p’-p  AeQ„  P  P  ~  P  P  “  To’To  "  P  P  " 

~  P 

(2.2-9) 


*See  Appendix  A. 


r’T«S' 


where  A  =  (£,m)1  is  a  vector  with  integer  components  and  the 

set  Q„  is  defined  as 
P 


Qp  =  {A|A^0,  | |s+epJp  A | |  <  2 | | s' | | } 


where 


s'  =  (  1/2t  j  ,  l/2ipT 


(2.2-10) 


(2.2-11) 


is  the  vector  of  Nyquist  frequencies  associated  with  the  grid 
M0  and  I  I*  I  I  denotes  vector  magnitude. 

p 

Equation  2.1-5  is  used  to  determine  the  statistics  of 
residuals  in  quantities  other  than  the  surface  anomalous  poten¬ 
tial.  If  w  is  a  vector  whose  components  are  arbitrary  field- 
related  quantities,  the  spectral  density  matrix  of  the  post¬ 
survey  residuals  in  w  is  obtained  from 


<t>  =  G  G  4> 

6w,<5w  -  -  6To,6Tq 


(2.2-12) 


where  G  is  the  vector  transfer  function  from  surface  anomalous 
potential  to  w.  The  inverse  Fourier  transform  of  <l>6w  yields 
the  covariance  matrix  of  the  residuals: 


00 

r6w.«w<£>  =  dsids: 


(2.2-13) 


from  which  root-mean-square  (rms)  values  of  the  residuals  in 
the  components  of  w  are  easily  obtained  by  taking  square  roots 
of  the  diagonal  entries  in  Rfiw  5w(0). 

In  many  instances  it  is  convenient  to  express  post¬ 
survey  residual  gravity  in  terms  of  spatially-averaged  values. 
Next,  the  formulas  for  the  average  spectral  density  of  the 


2-11 


1  -  '  »  *  _  »  .  »  * 
•  '«  •  k  *  h 


<*  •  »  •  »  '  J  «  * 4 


WWW 


residuals  are  particularized  to  the  estimation  of  spatially- 
averaged  quantities. 


Consider  a  square  block  with  sides  of  length,  a,  paral 

T 

lei  to  the  east  and  north  directions.  Let  x  =  (x^,  x2)  be 
the  center  of  the  block.  The  a-mean,  y  ,  of  a  field-related 
quantity  y  at  the  point  x  is  the  average  of  y  on  this  block, 
i  .  e .  , 

x2+a/2  _x1+a/2 

ya(x)  =  I  j  y(x')  dx|dx£  (2.2-14) 

a  x2-a/2  x^-a/2 

T 

where  x'  =  (xj ,  x^)  • 

Now,  suppose  that  the  vector  of  quantities  to  be  esti¬ 
mated  from  the  survey  data  consists  exclusively  of  a-raeans. 
Let  wa  be  the  field  vector  to  be  estimated.  It  is  shown  in 
Appendix  C  that  the  vector  transfer  function  from  anomalous 
surface  potential  to  the  estimated  vector  is  h  G  where  G  is 
the  vector  transfer  function  from  anomalous  surface  potential 
to  the  vector  w  whose  a-means  are  w  and 


h(s) 


sinnas^  sinnas2 


(2.2-15) 


is  the  transfer  function  associated  with  the  a-mean  averaging 
operation.  Consequently,  from  Eq.  2.1-5,  the  average  spectral 
density  of  the  residuals  in  the  a-means,  ^g^a  g^a,  is 


^Sw^Sw3  -  -  *61  ,6T 

—  —  o  o 


(2.2-16) 


with  g,j.  given  by  Eq.  2.2-7 
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From  Eqs.  2.2-12  and  2.2-16,  the  spectral  density 
matrix  of  a-mean  residuals  can  be  expressed  as 


4> 


6wa 


6w ,  6w 


(2.2-17) 


Thus,  to  evaluate  the  statistics  of  the  residuals  in  the  a-means 
it  suffices  to  multiply  the  spectral  density  of  the  residuals  in 
(the  point  quantities)  w  by  the  squared  magnitude  of  the  averag- 

A 

ing  transfer  function  h  as  in  Eq.  2.2-17. 


2.3  SURVEY  GEOMETRY  AND  MEASUREMENT  ERROR  MODELS 

This  section  presents  the  geometry  and  the  error  models 
associated  with  the  following  types  of  survey: 

•  Satellite  Radar  Altimetry 

•  High- Low  SST 

•  Land-based/Shipborne  Gravimetry 

•  Airborne  Gradiometry 

in  subsections  2.3.1  through  2.3'. 4,  respectively.  Survey  geom¬ 
etry  refers  to  the  position  of  the  measurements  and  the  orien¬ 
tation  of  the  sensors  with  respect  to  the  earth.  Measurement 
error  models  refer  to  the  enumeration  and  characterization  of 
the  various  error  sources  which  affect  the  individual 

measurements . 

The  geometry  and  type  of  survey  determine  the  rotation 

and  spacing  matrices  0„  and  J  and  the  transfer  functions  F„ 

P  P  P 

from  anomalous  surface  potential  to  the  quantities  being  meas¬ 
ured.  These,  in  turn,  characterize  completely  the  aliasing 
errors  once  a  field  model  is  specified  (Eq.  2.2-9). 
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The  measurement  error  models  determine  the  spectral 

densities  j,  of  Eq.  2.2-8.  Explicit  formulas  for  these 

-P’"P 

spectral  densities  are  given  in  this  section. 

2.3.1  Satellite  Radar-Altimeter  Survey 

The  geometry  of  a  satellite  radar-altimeter  survey  is 
illustrated  in  Fig.  2.3-1.  Satellite  groundtracks  are  viewed 
as  two  collections  of  parallel  equally- spaced  straight-line 
tracks  containing  data  at  regular  intervals  along  each  track. 

Two  measurement  classes  are  associated  with  the  survey: 
one  class,  A,  corresponds  to  the  ascending  tracks;  the  other,  D, 
corresponds  to  the  descending  tracks.  Let  the  grids  associated 
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DESCENDING  ASCENDING 

GROUND  TRACKS  GROUND  TRACKS 

r  A  i  r  A  "  — V 


Figure  2.3-1  Satellite  Altimeter  Survey  Geometry 
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with  these  two  classes  be  and  M^.  Contiguous  parallel 
tracks  are  separated  by  a  distance  x  in  the  east  direction 
and  successive  measurements  along  a  track  are  separated  by  a 
distance 


r. 


Gm. 


-.1/2 


Ta  ta  L(R+h)3. 


(2.3-1) 


where 
(t  - 


t  is  the  time  interval  between  altimeter  samples 
0.1  sec  for  GEOS-3  and  SEASAT-1),  Gmg  is  the  product  of 


the  gravitational  constant  and  the  mass  of  the  earth 


(Gm, 


3.986xl014  m3/sec2 ) , 


R  is  the  radius  of  the  earth 


(R  -  6.378x10°  m)  and  h  is  the  height  of  the  satellite  over 
the  earth's  surface  (h  -  8.0xl03  m)  for  GEOS-3  and  SEASAT-1. 

The  angle  a  is  the  angle  the  ascending  groundtracks 
form  with  the  east  direction  at  the  center  of  the  region  where 
estimates  are  sought  (the  origin  of  the  plane  of  the  flat-earth 
approximation).  This  angle  can  be  computed  from 


tan  a  = 


(cos2® 


cosy ) 


1/2 


cos  y  -  0.058834  cos*®  (^) 


R+h,  3/2 


(2.3-2) 


where  ®  is  the  latitude  of  the  origin  of  the  plane  of  the  flat- 
earth  approximation  and  y  is  the  inclination  of  the  satellite's 
orbit  (y  =  115°  for  GEOS-3  and  y  =  108°  for  SEASAT-1). 

The  grids  and  have  axes  parallel  to  those  of 
primed  reference  frames  (see  Fig.  2.2-1)  which  correspond  to 
rotation  angles  of  a  and  -a  with  respect  to  the  east  direction. 
Let  ©A  and  0^  be  the  rotation  matrices  of  the  grids  MA  and 
respectively.  These  matrices  are  given  by 


and 


(cos  or 
sin  a 


-sin  a 
cos  a 


(2.3-3) 


(cos  a 
-sin  a 


sin  a 
cos  a 


(2.3-4) 


The  spacing  matrices  JA  and  of  the  grids  MA  and 
coincide  with  each  other  and  are  given  by 


where  z  is  given  by  Eq.  2.3-1  and 

a 

x  =  t  sin  a 
c  e 


(2.3-5) 


(2.3-6) 


The  value  of  (Fig.  2.3-1)  can  be  found  from  the  equatorial 
separation  of  the  groundtracks ,  t  ,  by  the  simple  formula 

xe  =  Teq  cos  8  (2.3-7) 

Satellite  radar-altimetry  data  furnish  values  of  the 
undulation  of  the  geoid  corrupted  by  measurement  noise.  Thus, 
the  transfer  functions  from  anomalous  surface  potential  to  the 
quantities  being  measured  in  classes  A  and  D,  FA  and  F^,  are 
identical  to  each  other  and  correspond  to  the  transfer  function 
from  anomalous  surface  potential  to  undulation  of  the  geoid 
given  in  Table  2.1-1;  i.e., 


Fa(s)  =  FD(s)  =  l/gQ 


(2.3-8) 
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Measurement  error  models  are  discussed  next.  For  sym¬ 
metry  reasons,  the  spectral  densities  of  the  errors  in  classes 
A  and  D  agree  with  each  other  when  they  are  expressed  in  terms 
of  the  primed  coordinate  systems  of  their  respective  classes. 
Thus,  an  expression  for  the  spectral  density  of  the  errors, 

<1>E  £(s'),  in  terms  of  frequencies  s'  =  (s^.s^,)1  measured  in 
the  along-track  (s  )  and  cross-track  (s  )  directions  is  derived 
first.  The  spectral  densities  of  the  errors  in  classes  A  and 
D  are  then  obtained  in  terms  of  east  and  north  frequencies  by 
suitable  rotations  of  coordinate  systems. 

The  Nyquist  frequencies  in  the  along-track  and  cross¬ 
track  directions  are  l/2t  and  l/2x  ,  respectively.  Since  the 
measurement  error  spectral  density  4>i  -(s')  is  periodic  in  the 

L  ,  III  — 

along-track  and  .cross-track  directions  with  a  period  of  twice 
the  Nyquist  frequency  in  each  direction,  it  suffices  to  speci¬ 
fy  -(s')  on  the  domain  -l/2t  <s  <l/2r  and  -1/2t  <s  <1/2t  . 

L  ^  L  333  CCC. 

The  error  E(Q)  in  the  measurement  at  the  point 

T 

Q  =  (j,k)  of  any  one  of  the  grids  or  M^  is  taken  to  be  the 
additive  combination  of  three  independent  error  sources: 

E(Q )  =  N(Q )  +  C(fi)  +  B(fi)  (2.3-9) 

The  terms  N,  C  and  B  correspond  respectively  to 

•  Instrument  noise  and  sea-state  effects 
(scattering  of  the  radar  pulse  by  ocean 
surface  waves) 

•  Uncorrected  ocean-current  dynamic  height 

•  Post-adjustment  bias-type  orbit  and  tide- 
correction  errors. 


The  spectral  density  of  the  measurement  errors,  is  the 
sum  of  the  spectral  densities  of  the  terms  on  the  right-hand 
side  of  Eq.  2.3-9;  i.e., 


d>  • 

E,E 


=  <t> ' 
VN,N 


<t>  * 

C,C 


d>  > 

*B,B 


(2.3-10) 


First,  consider  the  instrument  noise  and  sea-state  er¬ 
rors  N.  These  errors  are  independent  from  measurement  point  to 
measurement  point  and  from  track  to  track.  Their  covariance  is 


RN,N^— ^  aN  ^  —  ^ 


(2.3-11) 


where  is  the  standard  deviation  of  the  error  N  in  each  meas¬ 
urement  (aN  =  0.6  m  for  GEOS- 3  and  0.1  m  for  SEASAT-1)  and 
where  6(0)  is  given  by 


l  1  if  0  =  0 

6(0)  =  j  (2.3-12) 

(  0  otherwise 


The  spectral  density  N  is  the  finite  Fourier  transform  of 
RN,N: 

,N(- '  }  =  e~i2n< ^a-’-' '  (2.3-13) 


where  JA  =  Jg  is  the  spacing  matrix  of  the  grid  and  A(JA)  =  x ax 
is  the  determinant  of  JA<  Since  jj(0)  vanishes  for  all  0/0, 
the  spectral  density  <t>^  ^  can  be  easily  evaluated.  The  result 
is  the  constant  (white)  spectrum 

*N,N(-' )  =  TaTcaN  (2.3-14) 

Next,  consider  the  uncorrected  ocean-current  dynamic 
height  errors  C(0).  Reference  4  presents  a  statistical  model 


2-18 


mvwxim  ^  ^  >  V’  V’*:“  T 


•J'm^  rr 


^  1 1 '  ■> '  K.i  v*  t« rrr^r  y  ^  »  «  7  ■■? 


for  the  spatial  distribution  of  the  ocean-current  induced  sea- 
surface  height  over  the  geoid  as  a  continuous  scalar  process 
defined  in  two-dimensional  space,  h(x').  The  spectral  density 
of  h  is  that  of  an  isotropic  third-order  Markov  model  (Ref.  5) 
given  by 


lOnOuPu 

d>  i  (s'  )  -  _ — _ - _ 

h-h  -  ’ 


(2.3-15) 


where  (in  m)  is  the  standard  deviation  of  the  sea-surface 
height  given  as  a  function  of  latitude  by 


=  0.677  sin  S  (2.3-16) 

and  where  1/0^  =  4.2*10  m  is  the  characteristic  distance  of 

the  model.  In  Eq.  2.3-15,  s’  =  Ms’ll  =  (s2+s2)1/2. 

Si  c 

The  spectral  density  4*^  ^  is  an  instantaneous  model 
in  the  sense  that  it  characterizes  the  spatial  variability  of 
the  sea-surface  height  at  a  fixed  time.  Actually  (Ref.  6), 
ocean-current  induced  sea-surface  height  does  not  vary  sub¬ 
stantially  for  time  spans  of  the  order  of  one  day.  However, 
after  a  period  of  one  or  two  weeks,  there  is  no  significant 
correlation  between  the  sea-surface  height  at  the  beginning 
and  end  of  the  period.  Consequently,  the  sea-surface  height 
model  of  Eq.  2.3-15  can  be  used  to  infer  the  behavior  of  the 
errors  C(Q)  for  a  single  track  of  data  but  cannot  be  used  for 
relating  the  errors  in  different  tracks. 


*The  characteristic  distance  1/0^  should  not  be  confused  with 
the  correlation  distance.  For  a  third-order  Markov  process 
the  correlation  distance  is  2.903/0^  (Ref.  5).  The  correlation 
distance  of  the  sea-surface  height  model  is  122  km. 
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In  practical  situations,  most  pairs  of  parallel  ground 
tracks  chosen  arbitrarily  over  the  region  where  estimates  are 
sought  correspond  to  satellite  passes  interspaced  by  periods 
of  several  weeks.  Thus,  the  ocean-current  induced  sea-surface 
height  errors  on  different  tracks  are  modeled  as  being  inde¬ 
pendent  of  one  another. 

The  spectral  density  of  the  sea-surface  height  on  a 

■k 

single  groundtrack  ,  ^(sa)  °btained  by  integrating  the 

two-dimensional  spectral  density  ^(s')  in  the  crosstrack 
direction  (see  Appendix  A,  Eq.  A-106): 


00 

Sh,h<sa>  s  f  4’h,h<5'>  ds, 


The  result  of  this  integration  is 


(2.3-17) 


Sh,h(sa^  ~ 


16o^3 

IPh+<2"sa>2I: 


(2.3-18) 


This  spectral  density  describes  the  behavior  of  the  groundtrack 
sea-surface  height  as  a  continuous  (one-dimensional)  process. 
The  spectral  density  of  the  measurement  errors  along  a  track 
of  data,  Sc  c(sa)  correspond  to  the  aliased  version  of  this 
continuous  process.  In  Appendix  C.2,  the  effects  of  sampling  a 
two-dimensional  process  are  analyzed.  A  similar  analysis  for 
one-dimensional  processes  yields  the  following  relation  between 
the  spectral  density,  ^  and  S^.  ^  of  the  continuous  and  sam¬ 
pled  versions  of  the  process 


*The  uppercase  letter  S  is  used  to  denote  one-dimensional 
spectral  densities. 


d.  .  V.  L.  -  *  .  V*.  1.  ' 


.  1 _ _ 
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(2.3-19) 


SC,C(sa)  =  2^  Sh,h  (sa+£/ta) 
£=-00 


Since  sea-surface  height  errors  in  different  tracks  are  inde¬ 


pendent,  the  (two-dimensional)  error  spectral  density  ^(s ' ) 


is  white  in  the  crosstrack  direction.  Thus, 


?C,C(2'>  =  'c  £  Sh.h<Vi/xa) 


(2.3-20) 


£  =  -oo 


A  graphical  interpretation  of  Eq.  2.3-19  is  presented 
in  Fig.  2.3-2.  The  PSD  of  the  errors  along  a  track  is  obtained 
as  the  infinite  superposition  of  translates  of  the  PSD  of  the 
continuous  process  representing  the  sea-surface  height.  Con¬ 


sider  the  interval  -l/2x  <s  <l/2x  .  Except  for  the  regions 

«  »  a 


near  the  edges  of  the  interval,  there  is  negligible  contribu- 

-t- 

tion  to  the  sum  in  Eq.  2.3-19  from  terms  for  which  I/O ' .  In 
the  regions  near  the  edges  of  the  interval,  the  PSD  of  the  er¬ 
rors  arising  from  instrument  noise  and  sea-state  effects, 
is  much  larger  (ten  orders  of  magnitude)  than  Thus,  for 

all  practical  purposes,  the  PSD  of  the  ocean-current  induced 


sea-surface  height  errors  in  the  interval  -l/2x  <s  <l/2x  can 

3  3  3 


be  taken  as 


C,CV- 


(s '  )  = 


16,00^5/3 


(pj  +  (27TSar] 


(2.3-21) 


For  values  of  s  outside  this  interval,  $1  r(s')  is  obtained 

3  L  j  L. 


by  repeating  Eq.  2.3-21  periodically  as  indicated  by  the  solid 
line  in  Fig.  2.3-2. 


*This  is  because  the  characteristic  distance  1/p^  is  two  orders 
of  magnitude  larger  than  the  sampling  spacing  x". 
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2-21 


«v  A.'i S-V ,V.'j  - .  -  .  -A.>  A 


Figure  2.3-2 


Graphical  Interpretation  of  Eq.  2.3-19 


Next,  consider  the  bias- type  orbit  and  tide  correction 
errors,  B(fl ) .  On  any  finite- length  of  single-track  data,  these 
errors  manifest  themselves  as  an  apparent  bias  with  a  standard 
deviation  of  ctr  =  0.5  m.  In  the  cross-track  direction,  the  er¬ 
rors  B(Q)  are  modeled  as  being  independent  from  track  to  track 
for  the  same  reasons  uncorrected  ocean-current  induced  sea- 
surface  height  errors  are  taken  as  independent  in  that  direction 

The  spectral  density  of  B(fi)  is  modeled  as  being 
Gaussian  in  the  along-track  direction  with  a  correlation  dis¬ 
tance  equal  to  the  radius  of  the  earth,  R,  and  as  being  white 
in  the  cross-track  direction.  Thus,  for  -l/2t  <s  <1/2t 

3  a  <3 

,  -*W 

tc  R°b  e  a  (2.3-22) 

For  |s  |  >  1/2t „  the  value  of  D(s  )  is  found  by  reproducing 

a  a  D ,D  3 

periodically  the  values  on  the  right-hand  side  of  Eq.  2.3-22. 

The  satellite  radar-altimeter  survey  error  spectrum 
*E  e^-'^  &iven  ky  the  sum  of  the  (periodically-continued) 


spectra  of  Eqs.  2.3-14,  2.3-21  and  2.3-22.  In  terms  of  fre¬ 
quencies  measured  in  the  east  (s1)  and  north  (s2>  directions, 
the  class  A  measurement  error  spectrum  can  be  computed  from 


,E^  — ^  ~ 


*E  , E (0A— ^ 


(2.3-23) 


and  the  measurement  error  spectrum  of  class  D  can  be  obtained 
from 

?E.E<S)  =  *i,E<8iU>  (2.3-24) 

where  ©A  and  ©D  are  given  by  Eqs.  2.3-3  and  2.3-4  and  where 
s  —  (s^  ,s2)  • 

2.3.2  High-Low  SST  Survey 

The  convention  on  the  position  of  the  measurements  of 
a  survey  consisting  of  high-low  SST  data  is  shown  in  Fig.  2.3-3 
The  survey  data  consist  of  range-rate  measurements  from  a  satel 
lite  in  synchronous  orbit  to  a  satellite  in  a  lower  orbit  cor¬ 
rected  for  the  nominal  motion  of  the  low  satellite  and  the 
spurious  motion  of  the  high  satellite.  In  other  words,  the 
data  consist  of  noisy  measurements  of  the  line-of-sight  com¬ 
ponent  of  the  variational  velocity  induced  on  the  low  satel¬ 
lite  by  the  gravity  disturbance  at  height  h. 

As  in  the  geometry  of  the  satellite  radar-altimeter 
survey,  measurements  on  the  earth  plane  lie  at  regular  inter¬ 
vals  on  equally-spaced  ascending  and  descending  groundtracks . 
The  same  notation  used  in  Subsection  2.3.1  is  adopted  here  for 
convenience.  Thus,  Eqs.  2.3-1  through  2.3-7  which  describe 
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Figure  2.3-3  Convention  on  the  Position  of 

High-Low  SST  Measurements 

the  geometry  of  the  grids  M.  (ascending  groundtracks )  and  Mn 

A  JL  U 

(descending  groundtracks)  hold  in  this  case  as  well.” 

The  transfer  functions  from  anomalous  surface  poten¬ 
tial  to  the  measurements  in  classes  A  and  D  are  considered 
next.  The  acceleration  perturbation,  v,  acting  on  the  low 
satellite  at  any  point  along  its  trajectory  is  the  gravity 
disturbance  vector;  i.e., 

v  =  (rXl ’rx2’r2)  (2.3-25) 


*In  Eq .  2.3-1,  t  becomes  the  time  interval  between 

a 

range-rate  samples  usually  taken  as  t  =10  sec. 


successive 


Therefore,  the  velocity  perturbation,  v,  is  the  time-integral 
of  the  gravity  disturbance  at  height  h  over  the  satellite's 
groundtrack. 


Let  V  be  the  speed  at  which  the  groundtrack  is  swept 


V  -  R 


<R+h>' 


(2.3-26 


Since  integration  with  respect  to  time  can  be  replaced  by  inte¬ 
gration  with  respect  to  distance  in  the  along-track  direction 
normalized  by  V,  the  vector  transfer  function  from  anomalous 
surface  potential  to  velocity  perturbation,  F  ,  is  the  product 
of  the  transfer  function  from  anomalous  surface  potential  to 
gravity  disturbance  at  height  h  and  the  transfer  function  asso¬ 
ciated  with  the  operation  of  spatial  integration  in  the  along- 
track  direction  normalized  by  the  constant  V.  It  is  shown  in 
Appendix  A  (Eq.  A-13)  that  the  transfer  function  corresponding 
to  integration  in  the  along-track  direction  is  l/(i2ns  )  where 
s  is  frequency  measured  in  the  along-track  direction.  Thus, 

cl 

from  Table  2.1-1 


-2nsh 


^v  = 


(s1 y  s  2 » is ) 


(2.3-27) 


At  any  point  along  the  satellite's  path,  the  line-of- 

sight  component  of  the  velocity  v  is  the  scalar  product  <u,v> 

where  u  is  a  unit  vector  in  the  direction  from  the  high  to  the 

low  satellite.  Measurement  geometry  is  illustrated  in 

Fig.  2.3-4.  The  high  satellite  lies  on  the  equatorial  plane 

at  a  distance  h  =  3.5786*10^  m  from  the  surface  of  the  earth, 
s 

The  low  satellite  is  at  latitude  9  and  at  longitude  A  measured 
with  respect  to  the  meridian  on  which  the  high  satellite  remains 
stationary.  The  unit  vector  in  the  direction  from  the  high  to 
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HIGH 

SATELLITE 


the  low  satellite  is  given  in  an  east-north-vertical 
the  position  of  the  low  satellite  by 


lu  l\ 

_  i 
p 

1  (R+hg)sinA  ^ 

u  = 

u2  1 

(R+h  )  sinScosA 
|  s 

^(R+h)-(R+hs)cos3cosA 

frame  at 


(2.3-28) 


where 

p  =  ( (R+h)2  +  (R+hs)2  -  2(R+h) (R+hg )cos$cosA ] 1/2  (2.3-29) 


T 

The  unit  vector,  u  =  (u^,U2,uz)  ,  varies  as  the  low 
satellite  changes  its  position  over  the  earth.  This  implies 
that  the  measurement  equations  are  time-varying.  However, 
gravity  estimates  are  sought  on  a  small  region  of  the  earth 
over  which  the  unit  vector  u  can  be  approximated  by  a  constant 
vector.  This  constant  vector  is  chosen  as  that  obtained  from 
Eq.  2.3-28  when  3  and  A  are  interpreted  as  the  coordinates  of 
the  origin  of  the  plane  of  the  flat-earth  approximation. 
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The  transfer  function  from  anomalous  surface  potential 
to  measured  quantities  is  given  by  the  scalar  product  <u,F  > . 
For  the  measurements  of  class  A,  the  along-track  frequency  s 

<3 

appearing  in  Eq.  2.3-27  can  be  expressed  in  terms  of  east  (s^) 
and  north  (S2)  frequencies  as 


=  s^  cosa  +  S2  sina  (2.3-30) 

where  a  is  the  angle  of  Eq .  2.3-2.  Consequently,  the  transfer 
function  of  class  A  is 


(u, s, +u0s0+iu  s)e 

p  (s\  -  — ± — ± ± — - _ ? _ 

A  —  V(s^cosa  +  S2sina) 


Similarly,  the  transfer  function  of  class  D  is 


(u,s,+u~s9+ius)e 

F  ( s )  -  — - — ± — = — - £ - 

D  — '  V(s^cosa  +  S2sina) 


(2.3-31) 


(2.3-32) 


A  simplified  measurement  error  model  has  been  used  to 
obtain  the  results  of  Chapter  4.  The  errors  in  the  measure¬ 
ments  are  modeled  as  being  independent  of  one  another.  The 
spectral  densities  of  the  errors  in  classes  A  and  D  are  iden¬ 
tical  to  each  other  and  are  given  by 


*E , E^-^  ~ 


x  t  a 
a  c 


(2.3-33) 


where  x  and  x  are  the  measurement  spacings  in  the  along-track 

3  C 

and  cross-track  directions  and  where  o  is  the  standard  devia¬ 
tion  of  the  error  in  each  measurement. 
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2.3.3  Land-Based/Shipborne  Gravimetric  Survey 

Land-based  and  shipborne  gravimetric  surveys  are  con¬ 
sidered  in  this  subsection.  Shipborne  gravimetric  survey  geom¬ 
etry  and  measurement  error  models  are  discussed  first.  Corre¬ 
sponding  models  for  land-based  surveys  are  then  obtained  by 
removing  from  the  shipborne  gravimetric  survey  models  those 
effects  which  are  particular  to  ocean  surveys. 

Shipborne  Gravimetric  Survey  -  The  geometry  of  a  ship¬ 
borne  gravimetric  survey  is  illustrated  in  Fig.  2.3-5.  Survey 
data  are  collected  at  regular  intervals  by  a  gravimeter  on  board 
a  ship  traveling  at  constant  speed  along  parallel  equally-spaced 
east-west  tracks.  Since  the  measurement  grid  is  oriented  with 
the  east-north  reference  frame,  the  rotation  matrix,  0,  is  the 
identity,  I.  The  spacing  matrix  is 


R-40096 


T 

r2 
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Figure  2.3-5 


Shipborne  Gravimetric  Survey 
Measurement  Geometry 


where  is  the  spacing  between  tracks  and  where 

t1  =  V  At  (2.3-35) 


is  the  spacing  in  the  east  direction  given  in  terms  of  the 
nominal  ship’s  speed  V  and  the  time  interval  between  samples 
At.  Typical  values  are  assumed  for  x^,  V,  and  At  based  on 
NAVOCEANO's  survey  patterns. 

The  transfer  function  from  surface  anomalous  poten¬ 
tial  to  the  quantity  being  measured  is  discussed  next.  Denote 
by  £  the  difference  between  the  gravimeter  readings  and  refer¬ 
ence  gravity  at  the  points  where  measurements  are  taken.  Be¬ 
cause  of  the  motion  of  the  ship  over  the  earth's  surface,  the 
quantity  £  is  not  the  gravity  anomaly.  To  account  for  the 
ship's  motion,  several  corrections  are  applied  to  the  acceler¬ 
ation  £  . 


T 

Let  v  =  (v^.v2^  and  ®  be  the  true  velocity  of  the 

survey  ship  with  respect  to  the  earth  and  its  true  latitude. 

The  gravity  field  includes  the  centrifugal  acceleration  due  to 

the  rotation  of  the  earth.  This  acceleration  is  directed  away 

2 

from  the  axis  of  rotation  and  has  a  magnitude  Ru>  cos&  at  lati- 

-4  c 

tude  S  (uj=0 . 7292115147*  10  rad/sec).  The  component  directed 

t  2  2 

towards  the  center  of  the  earth  is  -Ru>  cos  St.  The  motion  of 
the  ship  in  the  east  direction  causes  the  instantaneous  angular 


*E(v1)  =  V,  E(v2)  =  0 


velocity  to  differ  from  the  angular  velocity  of  the  earth  by 
an  amount  6ui  given  by 


Suj 


V1 

R  cos8t 


(2.3-36) 


To  first-order,  the  corresponding  contribution  to  the  gravimeter 
measurement  is 


=  -2Ruj  cos^d^  6u)  (2.3-37) 

or,  combining  Eqs .  2.3-36  and  2.3-37, 

t;  '  =  -2u»v^  cos&t  (2.3-38) 

In  addition  to  the  contribution  resulting  from  the  difference 
in  rotation  rates,  there  is  an  acceleration  of 

V'  =  -<vj  +  v|)/R  (2.3-39) 

resulting  from  the  motion  of  the  ship  over  a  curved  surface. 

If  the  gravimeter  readings  were  error-free,  the  gravity 
anomaly  could  be  computed  from 

Ag  =  £  -  C  -  C"  (2.3-40) 

However,  the  velocity  and  latitude  which  enter  into  the  compu¬ 
tations  in  Eqs.  2.3-38  and  2.3-39  are  not  perfectly  known.  The 
evaluation  of  Eq.  2.3-40  is  performed  using  the  Inertial  Navi¬ 
gation  System  ( INS )- indicated  velocity  v°  =  (v°,  v^)^"  and  lati- 
tude  d°.  Thus,  shipborne  gravimetric  data  consist  of  the 


following  combination  of  gravimeter  readings  and  measurements 
from  an  inertial  navigator: 

4>  =  £  +  2wv°cos3°  +  ||v°||2/R  (2.3-41) 

By  adding  and  subtracting  the  true  correction  £'  +  £" ,  Eq .  2.3-41 
can  be  written  as 

4»  =  (£-£'-£")  +  2iucos3°  6v^  +  2u»v°sin3°  63 

+  2v°6v1/R  +  2v°6v2/R  (2.3-42) 

where  second-order  terms  have  been  neglected  and  where 

6v  =  (6v1,6v2)T  =  (v^-Vj^,  v2-v2 ^  (2.3-43) 

is  the  error  in  the  velocity  estimate  from  the  INS  and 
63  =  3°  -  3 ^  is  the  error  in  latitude. 

When  typical  values  of  6v^  ,  6v2 ,  and  63  are  considered, 
the  combined  contribution  of  the  last  three  terms  in  the  right- 
hand  side  of  Eq.  2.3-42  is  several  orders  of  magnitude  smaller 
than  that  of  the  term  2u>cos3°  6v^ .  Therefore,  the  last  three 
terms  in  Eq.  2.3-42  can  be  neglected.  On  the  other  hand,  since 
the  latitude  3°  does  not  vary  substantially  over  the  survey 
region,  the  angle  3°  can  be  replaced  by  3,  the  latitude  of  the 
origin  of  the  plane  of  the  flat-earth  approximation.  Thus, 
for  all  practical  purposes,  Eq.  2.3-42  can  be  rewritten  as 

=  £-£•-£"  +  2u>cos3  6v^  (2.3-44) 

The  term  2wcos3  6v^  is  called  the  Eotvos  correction 
error.  It  is  caused  by  the  error  in  the  estimate  of  the  east 
component  of  velocity  obtained  from  the  INS  aboard  the  survey 
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ship.  The  east-velocity  error,  in  turn,  arises  from  the  east 
component  of  the  gravity  disturbance  vector  sensed  by  the  iner¬ 
tial  navigator  and  from  random  accelerometer  and  gyro  errors. 
Consequently,  if  shipborne  gravimetric  data  are  viewed  as  grav¬ 
ity  anomaly  measurements,  the  measurement  errors  are  correlated 
with  the  gravity  field.  This  correlation  must  be  accounted 
for  in  data  processing  algorithms. 

An  equivalent  formulation  is  obtained  when  the  data 
are  viewed  as  consisting  of  noisy  measurements  of  linear  combi¬ 
nations  of  the  gravity  anomaly  and  the  east  component  of  the 
gravity  disturbance  vector  modulated  by  the  INS  response. 

This  convention  is  adopted  here.  Let  u  and  y  be  the  INS  east- 
velocity  errors  induced  by  the  gravity  field  and  by  random 
errors  in  the  inertial  navigator  instruments  (accelerometer 
and  gyro  errors),  respectively.  Equation  2.3-44  becomes 

4*  =  (£  -  t '  -  £" )  +  (2wcos8)u  +  (2u>cos3)y  (2.3-45) 

The  term  (2uicosfl)u  is  considered  part  of  the  measured  quantities. 
Measurement  errors  include  the  term  (2iucosS)y.  In  this  formu¬ 
lation,  measurement  errors  turn  out  to  be  independent  of  the 
gravity  field. 

In  order  to  determine  the  transfer  function,  F,  from 
the  anomalous  surface  potential  to  the  quantity  being  meas¬ 
ured,  it  is  necessary  to  relate  the  velocity  errors  u  and  the 

east  component  of  the  gravity  disturbance  vector  r  .  The  INS 

X1 

transfer  function  from  sensed  acceleration  to  indicated  veloc¬ 
ity  is  very  well  approximated  by  that  of  a  damped  second-order 
system.  When  the  input  acceleration  and  the  output  velocity 
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are  expressed  as  functions  of  time  the  (steady-state)  transfer 

JU 

function” ,  G(f),  is  given  by 

G(f)  =  - (2.3-46) 

2pfof+i(f2-f2) 

where  f  is  frequency  measured  in  Hertz,  p  is  the  INS  damping 

1/2 

coefficient  and  f  =  (gQ/R)  '  /2n  is  the  Schuler  frequency 

(f  =  1.9708xl0‘4  Hz), 
o 

The  INS  transfer  function  from  sensed  acceleration  to 
indicated  velocity  as  functions  of  distance  in  the  east  di¬ 
rection,  G(s^),  can  be  obtained  by  making  the  transformation 
f  =  Vs^.  Thus,  G(s^)  =  G(Vs^),  and  from  Eq .  2.3-46, 

s1/(2nV) 

G( s,  )  =  - ± (2.3-47) 

1  2P(fo/V)Sl+i[sJ-(foA)2] 

It  then  follows  from  Eqs.  2.3-40  and  2.3-45  that  the  shipborne 
gravimetric  survey  transfer  function,  F,  is  given  by 

F(s)  =  [2ns-2/R]  +  2tucos$  [  i2ns^  ]G(  s^ )  (2.3-48) 

where  the  two  quantities  in  brackets  are  the  transfer  functions 
from  anomalous  surface  potential  to  gravity  anomaly  and  to  the 
east  component  of  the  gravity  disturbance  vector  at  the  surface, 
respectively . 

Note  that  the  relation  f  =  Vs^  used  in  deriving 
Eq.  2.3-47  is  an  approximation  because  the  east  velocity  is 


*The  integration  in  the  definition  of  Fourier  transform  is 
carried  out  over  the  time  domain. 
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not  truly  constant.  However,  its  use  in  obtaining  Eq.  2.3-47 
is  equivalent  to  neglecting  second-order  terms. 


The  shipborne  survey  measurement  error  model  is  dis¬ 
cussed  next.  Consider  first  the  uncorrelated  Eotvos  correction 
error,  C,  corresponding  to  the  term  (2uicosd)y  in  Eq .  2.3-45. 
These  errors  are  modeled  as  arising  from  a  white  noise  sensed 
acceleration  input  with  a  constant  spectral  density  q.  There¬ 
fore,  the  spectral  density  of  the  velocity  errors,  y ,  on  a 
single  track  of  data  is 


sy  ,y<sl>  =  q|G(s1) | 


(2.3-49) 


The  value  of  q  was  chosen  by  assuming  that  the  veloc¬ 
ity  errors,  y,  have  an  rms  value  oy  for  a  typical  INS,  i.e., 


2  _ 

CTy 


S  ( s, )  ds, 
y  ,y  1  1 


(2.3-50) 


This  results  in  a  value  of  q  =  8nVf  p a  . 

o  y 

Velocity  errors  arising  from  accelerometer  and  gyro 
errors  are  modeled  as  independent  on  different  tracks.  Fol¬ 
lowing  the  same  reasoning  used  in  obtaining  Eq .  2.3-21  from 
Eq .  2.3-18,  the  spectral  density  of  the  uncorrelated  Eotvos 
correction  errors  in  the  interval  -l/2t s^< l/2t ^  can  be  writ¬ 
ten  as 


C,C'- 


(s)  =  4uj^cos^8x  2<i  I  G(  s^  )  |  ^ 


(2.3-51) 


where  x^  and  X2  are  defined  in  Fig.  2.3-5 
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The  gravimeter  recordings  themselves  also  contain 

various  errors.  First ,  there  are  quantization ,  instrument 

noise  and  improperly  filtered  heave  motion  induced  errors. 

These  are  all  combined  into  a  single  error  term,  I,  modeled  as 

being  independent  from  point  to  point  and  track  to  track,  and 

“6  2 

having  an  rms  value  of  =  1.0x10  m/sec^  (0.1  mgal).  The 
PSD  of  I  is 

*I,I<s)  =  (2.3-52) 

Second ,  the  platform  is  not  directly  on  the  geoid  because  of 
tidal  and  ocean  current  dynamic  sea-surface  height  effects.  If 
the  platform  is  at  height  h  over  the  geoid,  the  gravimeter  read¬ 
ing  is  in  error  by  an  amount  -2gQh/R  (2gQ/R  =  0.308  mgal/m). 

In  the  open  ocean,  the  effects  on  the  gravimeter  readings  caused 
by  tidal  effects  can  be  predicted  to  an  accuracy  of  better  than 
0.03  mgal  (Ref.  7).  Consequently,  these  effects  are  neglected 
in  the  error  model. 


Consider  the  errors  K  induced  in  the  gravimeter  read¬ 
ings  by  ocean  currents.  Ocean  current  induced  sea-surface 
height,  h,  is  modeled  as  in  Eq.  2.3-15.  The  spectral  density 
of  h  (seen  as  a  continuous  process)  is 


Vh^)  = 


10"°h»h 


<2ns)2P/i 


(2.3-53) 


As  opposed  to  the  altimeter  survey  case,  successive  tracks  are 
surveyed  within  relatively  short  periods  of  time.  Thus, 

Eq.  2.3-53  is  adopted  directly  as  the  model  for  sea-surface 
height  in  the  survey  region.  The  spectral  density  of  the 
aliased  version  of  h,  H,  is  given  by 


<  s  +  .J-iA  ) 


(2.3-54) 


^  < 
11 ,11 


) 

A 


h  ,h 


where  the  summation  is  taken  over  all  integer-valued  vectors 
T 

A  =  ( j  ,k)  and  where  J  is  as  in  Eq .  2.3-34. 


For  frequencies  s  such  that  -l/2x  s^<  l/'2x  ^  and 

-l/2t2<S2<l/2t2  the  contribution  to  the  sum  in  Eq.  2.3-54  from 
all  terms  for  which  j/0  is  negligible  because  the  along-track 
Nyquist  wavelength,  2x^,  is  two  orders  of  magnitude  smaller 
than  the  correlation  distance  of  the  sea-surface  height  model 
(122  km).  However,  in  some  instances,  the  Nyquist  wavelength 
in  the  crosstrack  direction,  2^.  is  comparable  to  the  cor¬ 
relation  distance  of  the  sea-surface  height  model.  For 
-l/2x^<s^<l/2t^  and  - l/2x ^ s2< l/2x 2 >  the  aliased  spectrum  of 
the  sea-surface  height  is  approximated  as 

1 

$H,H(^)  =  S  *h,h  [s+(0,k/x2)T]  (2.3-55) 

k=-l 

For  the  same  range  of  frequencies,  the  corresponding  PSD  of 
the  gravimeter  measurement  errors,  K,  induced  by  sea-surface 
height  is 


®k,k(5> 


(2.3-56) 


with  ^  as  in  Eq .  2.3-55. 

In  the  absence  of  any  other  error  sources,  the  spectral 
density  of  the  shipborne  survey  measurement  errors,  r  is 

CL  9  CL 

given  by  the  sum  of  the  spectral  densities  of  the  velocity- 
induced  errors  (Eq.  2.3-51),  the  quantization  errors  (Eq.  2.3-52) 
and  the  ocean  current- induced  errors  (Eq.  2.3-56);  i.e., 


,  E(  — ^  =  4>C,C(-)  +  *1,1*-*  +  *K,K(-} 


(2.3-57) 


for  -l/2x s^< 1/2t ^  and  -l/2t2<S2<l/2t2 •  However,  actual  sur¬ 
veys  conducted  by  NAVOCEANO  follow  the  geometry  of  Fig.  2.3-5 
on  nearly  rectangular  blocks.  An  example  of  the  characteristic 
block  layout  is  presented  in  Fig.  2.3-6.  Often,  two  contiguous 
blocks  share  a  common  region  but,  in  some  instances,  there  are 
small  data  gaps  (islands,  reefs,  etc.)  in  the  coverage.  Since 
track  spacing  is  chosen  on  the  basis  of  high-frequency  energy 
content  in  the  gravity  field,  neighboring  blocks  tend  to  have 
the  same  track  spacing  but  the  tracks  in  one  block  are  not,  in 
general  the  continuation  of  tracks  in  another  block. 


R-62443 


Figure  2.3-6 


Characteristic  Block  Layout  of 
NAVOCEANO  Surveys 
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A  ronservat ive  approach  was  used  to  take  into  account 
these  irregularities:  Only  the  information  contained  in  ship- 
borne  survey  data  in  the  shorter  wavelengths  is  used  to  arrive 
at  gravity  estimates.  This  was  modeled  by  setting  the  measure¬ 
ment  error  spectra  4>^.  ^  equal  to  infinity  for  longer  wave¬ 
lengths.  Thus,  Eq.  2.3-57  was  modified  as 


(  r(s)  +  $,  (s)  +  $  (s) 

♦E  e<s>  =  ’  ’  ’  ” 

(  0O 


i  f  s  >  s 

o 

if  s  <  s 

o 


(2.3-58) 


2  2  1/2 

with  s- |  | s  |  | = ( s, +s9  )  and  s  equal  to  the  frequency  cutoff 

^  °  30 

selected.  For  numerical  computations  the  number  10  was  used 
in  place  of  infinity.  This  number  is  at  least  10  orders  of 
magnitude  larger  than  any  other  quantity  appearing  in  the  compu¬ 
tations  when  MRS  units  are  used. 


Summarizing,  for  a  shipborne  gravimetric  survey  the 
transfer  function  from  anomalous  surface  potential  to  the  quan¬ 
tities  being  measured  is  given  by  Eq .  2.3-48  and  the  spectral 
density  of  the  measurement  errors  is  given  by  Eq.  2.3-58. 


Land-Based  Gravimetric  Survey  -  The  model  for  the 
geometry  of  a  land-based  gravimetric  survey  is  illustrated  in 
Fig.  2.3-7.  Measurement  points  are  assumed  to  form  a  grid 
oriented  with  the  east  and  north  axes.  The  distances  between 


data  points  in  the  east  and  north  directions  are  t ^  and 
respectively.  It  then  follows  that,  as  in  the  shipborne  gravi¬ 
metric  survey,  the  rotation  matrix,  0,  of  the  measurement  grid 
is  the  identity  and  the  spacing  matrix,  J,  is  given  by  an  ex¬ 
pression  identical  to  Eq.  2.3-34.  Land-based  gravimetric  data 
are  thus  assumed  to  have  been  gridded  prior  to  processing. 
Research  is  currently  in  progress  to  account  directly  for  the 
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Figure  2.3-7  Land-Based  Gravimetric  Survey  Geometry 

irregular  geometry  associated  with  the  collection  of  land-based 
gravimetric  data. 

Data  consist  of  measurements  of  gravity  anomaly  re¬ 
duced  to  the  geoid.  The  transfer  function  from  anomalous  sur¬ 
face  potential  to  measured  quantities  is 

F(s)  =  2ns  -  2/R  (2.3-59) 

Measurement  errors  are  modeled  as  being  independent 
from  point  to  point  with  standard  deviation  a.  Their  spectral 
density  is 


~  2 
*^E  ,  E ^  ^  -  T  iT2a 


(2.3-60) 


Instrument  errors,  gridding  errors,  and  errors  in  the  process 
of  reducing  the  measurements  to  the  geoid  are  included  in 


Eq.  2.3-60.  For  the  results  given  in  Chapter  3  of  this  re¬ 
port,  Ti  =  T2  =  1-5*10^  m  and  a  =  3x10  ^  m/sec^  (3  ragal). 


2.3.4  Airborne  Gradiometer  Surve\ 


Figure  2.3-8  presents  the  convention  on  the  position 
of  the  measurements  of  an  airborne  gradiometer  survey.  Data 
are  collected  by  a  triad  of  gradiometers  aboard  an  aircraft 
flying  on  parallel  equally  spaced  east-west  tracks  at  an  alti¬ 
tude  h  above  the  surface  of  the  earth.  The  spacing  in  the 
north  direction  is  X2  and  the  spacing  in  the  east  direction  is 

tx  =  Vt  (2.3-61) 

where  V  is  the  (nominal)  speed  of  the  aircraft  and  t  is  the 
time  interval  between  successive  samples.  Thus,  the  rotation 
matrix  of  the  measurement  grid,  ©,  is  the  identity  and  the 
spacing  matrix  is 


Figure  2.3-8 


Airborne  Gradiometer  Survey 
Mecsurement  Geometry 


Prototype  gravity  gradiometers  developed  by  the  Bell 
Aerospace  Division  of  Textron,  Inc.  (Bell)  and  the  Charles 
Stark  Draper  Laboratory  (CSDL)  were  considered  in  the  analysis. 
Measurement  equations  for  airborne  surveys  using  the  Bell  and 
CSDL  gradiometers  are  given  below. 

Bell  Gradiometer  Survey  Transfer  Function  -  A  schematic 
diagram  of  a  Bell  gradiometer  is  presented  in  Fig.  2.3-9.  The 
instrument  uses  four  matched-accelerometers  mounted  on  a  slowly 
rotating  (0.25  Hz)  table.  The  outputs  of  the  accelerometers  are 
mixed,  preamplified,  band-limited  and  demodulated  at  twice  the 
rotation  frequency  to  yield  two  gradient  measurements,  one  at 
0  deg  phase  and  the  other  at  90  deg.  Let  the  axes  x' ,  y'  and 
z'  be  defined  as  in  Fig.  2.3-9.  One  of  the  measurements  (inline 
channel),  v^,  consists  of  the  difference  between  the  two  inline 
gradients  lying  on  the  plane  of  rotation  divided  by  two: 

v.  =  (T  ,  ,  -  T  ,  , )/2  (2.3-63) 

i  x'x'  yy 

R— 28144b 
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Figure  2.3-9 


Bell  Gradiometer  Idealization 


The  other  measurement  (cross-channel),  v  ,  yields  the  rotation 
plane  cross-gradient: 


vc  ~  *~x  '  y  ' 


(2.3-64) 


The  geometry  for  the  configuration  of  the  gradiometer 
triad  was  chosen  so  that  the  instruments'  spin  axes  coincided 
with  the  vertical  (z),  north  (X2)  and  east  (x-^)  directions 
(Fig.  2.3-10).  These  orientations  are  consistent  with  test 
data  provided  by  Bell  for  their  Baseline  gradiometer  instru¬ 
ment.  Consequently,  from  Eqs .  2.3-63  and  2.3-64  it  follows 
that  the  measurement  equations  are  given  by 


xlxl 


X2X2 


xlx2 


(2.3-65) 

where  v^,  v^  and  v^  are  the  inline-channel  outputs  and  V2 ,  v^ 
and  v6  the  cross-channel  outputs  from  the  gradiometers  whose 
spin  axes  are  oriented  in  the  vertical,  north  and  east  direc¬ 
tions,  respectively.  The  vector  transfer  function,  F,  from 
anomalous  surface  potential  to  the  measured  quantities  is  easily 
obtained  from  Table  2.1-1.  The  result  is 
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CSDL  Gradiometer  Survey  Transfer  Function  -  The  CSDL 
gradiometer  (Fig.  2.3-11)  employs  a  floated  electrostatically 
suspended  sphere  with  dense  material  packed  symmetrically 
about  two  opposite  poles.  Gravity  gradients  induce  a  torque 
on  the  sphere  that  is  sensed  through  the  voltages  needed  to 
hold  the  sphere  in  its  nominal  attitude.  Two  measurements  are 
obtained  from  a  single  gradiometer  instrument.  With  the  co¬ 
ordinates  x',  y'  and  z'  defined  as  in  Fig.  2.3-11,  the  meas¬ 
urements  are 


x ' 


=  r 


X  '  z 


V 


=  r 


y'z' 


(2.3-68) 
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Figure  2.3-11 


CSDL  Gradiometer  Float  Element 


The  geometry  considered  for  the  CSDL  gradiometer  triad 
is  shown  in  Fig.  2.3-12.  The  axes  of  weights  of  the  three 
gradiometers  (Zp  Z2  and  z^)  are  all  oriented  at  an  angle  of 
35.2644  deg  with  respect  to  the  vertical.  The  projection  of 
the  z^  axis  on  the  horizontal  plane  coincides  with  the  east 
direction  and  forms  angles  of  120  deg  with  the  projections  of 
the  axes  of  weights  of  the  other  two  gradiometers  on  the  same 
plane.  Let  and  V2  be  the  measurements  produced  by  the  gra¬ 
diometer  whose  weight  axis  is  z^ ,  and  those  produced  by 
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Figure  2.3-12 


CSDL  Gradiometer  Triad  (Tetrahedron 
Geometry) 


the  grad  iome  t  er  with  weight  axis  2.^ ,  and  v^  and  those  pro¬ 
duced  by  the  gradiometer  with  weight  axis  z^-  The  measurements 
can  be  expressed  in  terms  of  the  gradients  of  the  disturbance 
potential  in  the  east-north-vertical  frame  (Ref.  8)  by 
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(2.3-69) 

It  then  follows  that  the  vector  transfer  function  from  anoma¬ 
lous  surface  potential  to  measured  quantities  has  the  same 
form  as  that  obtained  for  the  Bell  triad;  i.e., 

F(s)  =  C  H  (s)  (2.3-70) 

where  C  is  the  6x6  matrix  appearing  in  Eq .  2.3-69  and  H(s)  is 
the  vector  transfer  function  from  anomalous  surface  potential 
to  the  gradients  of  the  anomalous  potential  at  height  h  given 
by  Eq.  2.3-67. 


Gradiometer  Survey  Measurement  Error  Models  -  The 
following  measurement  error  sources  were  considered  in  the 
simulations  performed: 


Instrument  noise  (self-noise) 


•  Errors  induced  by  the  mechanical  vibra¬ 
tions  to  which  the  gradiometers  are  sub¬ 
jected  aboard  an  aircraft. 

Spectral  densities  for  the  time  variation  of  the  self-noise  of 
the  Bell  Baseline,  Bell  Ball-Bearing  and  CSDL  gradiometers 
were  derived  from  the  manufacturers'  test  data  in  Ref.  9.  The 
same  self-noise  spectra  was  obtained  for  both  channels  of  each 
of  the  instruments  considered.  In  addition,  no  significant 
correlation  was  observed  between  the  two  outputs  of  each 
instrument . 

As  an  example  of  the  self-noise  spectra  derived  in 
Ref.  9,  Fig.  2.3-13  presents  the  spectra  of  the  self-noise  of 
the  Bell  Baseline  gradiometer  with  its  spin  axis  in  the  verti¬ 
cal  and  horizontal  positions .  The  shape  of  the  spectra  shown 
in  Fig;  2.3-13  is  typical  of  that  obtained  for  the  Bell  Ball- 
Bearing  and  CSDL  instruments  as  well.  The  spectrum  decays  at 
a  rate  of  6  db/octave  (red  noise)  and  flattens  out  (white  noise) 
at  high  frequencies.  The  analytic  form  of  such  a  spectral 
density  is 

Sr  jj(f)  =  +w  (2.3-71) 

where  f  is  measured  in  Hz.  For  the  Bell  Baseline  gradiometer 
the  values  of  r  and  w  are 

r  =  2.0xl0"6  E2  •  Hz 

w  =  81  E2/Hz 


(2.3-72) 
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Figure  2.3-13  Self-Noise  Spectra  of  the  Bell 

Baseline  Gradiometer 


when  the  spin  axis  is  vertical,  and 


r  =  16xl0'6  E2  •  Hz 


(2.3-73) 


w  =  86  E^/Hz 


when  the  spin  axis  is  horizontal  (Ref.  9).  For  the  Bell  Ball- 
Bearing  gradiometer,  data  corresponding  to  an  inclination  of 
the  spin  axis  with  respect  to  the  vertical  of  55  deg  were 
analyzed.  The  results  are 


r  =  7.7xl0"6  E2  •  Hz 


(2.3-74) 


w  =  290  E^/Hz 


These  values  were  used  for  both  horizontal  and  vertical  spin 
axes  of  the  gradiometers .  For  the  CSDL  gradiometers ,  the  values 
of  r  and  w  were  found  to  be 


r  =  0.2*10"6  E2  •  Hz 
w  =  2.3  E2/Hz 


(2.3-75) 


The  effects  of  translational  and  angular  vibrations  on 
the  gradiometer  measurements  are  also  analyzed  in  Ref.  9.  The 
net  result  is  an  increase  in  the  power  level  of  the  white  noise 
floor,  w,  which  depends  on  the  orientation  of  the  instrument. 
Vibration  data  corresponding  to  the  Bell  Baseline  instrument 
were  received  from  the  manufacturer.  The  sensitivities  obtained 
were  used  in  determining  the  modified  white  noise  levels  for  the 
Bell  instruments.  In  the  absence  of  vibration  data  from  CSDL, 
no  vibrationally  induced  errors  were  included  in  the  CSDL  gra¬ 
diometer  error  model. 


Because  of  the  way  the  instruments  sense  and  process 
the  gradients  of  the  field,  there  is  no  correlation  between 
vibrationally  induced  errors  in  both  channels  of  a  single  gradi¬ 
ometer.  Vibrationally  induced  errors  at  the  outputs  of  differ¬ 
ent  instruments  were  modeled  as  uncorrelated. 


Other  sources  of  error  which  were  not  included  in  the 
analysis  because  their  magnitude  is  very  small  compared  to  the 
self-noise  and  vibrationally  induced  errors  are: 


•  Thermally  induced  errors.  Temperature 

sensitivity  data  for  the  Bell  gradiometers 
(Refs.  10,  11  and  12)  indicate  that,  other 
than  the  red  noise  discussed  in  the  fore¬ 
going,  these  errors  can  be  neglected. 


/ 


•  Motion  correction  errors.  Airborne  gradi 
ometer  measurements  must  be  corrected  to 
account  for  the  fact  that  the  instruments 
are  not  stationary  in  inertial  space. 

The  corrections  are  on  the  order  of 

4  Eotvos.  Errors  in  the  corrections 
amount  to  a  few  hundredths  of  an  Eotvos. 

•  Errors  due  to  uncertainty  in  the  air¬ 
craft's  altimeter  readings.  These  er¬ 
rors  are  of  the  same  order  of  magnitude 
as  the  motion  correction  errors. 


In  addition,  the  gradiometers  are  naturally  susceptible  to 
time-varying  gravity  gradient  fields  caused  by  relative  motions 
of  nearby  masses  such  as  gimbals  in  the  inertial  platform. 
Corrections  can  be  applied  to  the  gradiometer  output  signals 
to  compensate  for  the  gradient  fields  caused  by  disturbing 
masses  whose  positions  can  be  monitored.  Errors  in  these  cor¬ 
rections  are  unaccounted  for  in  the  results  given  in  Chapter  3. 


The  spectral  density  of  the  survey  measurement  errors 
was  taken  to  coincide  on  a  single  data  track  with  the  spectral 
density  given  by  Eq .  2.3-71  (with  suitable  choices  for  the 
parameters  r  and  w)  after  transforming  time  frequencies  into 
spatial  frequencies  through  the  mapping 


f  =  V  sx  (2.3-76) 

where  V  is  the  aircraft's  speed.  On  different  data  tracks, 
measurement  errors  were  taken  to  be  independent.  The  measure¬ 
ment  error  spectral  density  matrix,  g(s)  ,  is  diagonal  with 
entries  along  the  diagonal  given  by 


E 


m 


(s) 


(2.3-77) 
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21 

(E  -Hz) 

2. OxlO'6 

7.7xlO'6 

2. OxlO"7 

r2 

(E2*Hz) 

2 .  OxlO"6 

7.7xl0‘6 

2. OxlO"7 

23 

(E  -Hz) 

16xl0“6 

7.7xl0‘6 

2. OxlO"7 

(E2^Hz) 

16xl0‘6 

7.7xlO-6 

2. OxlO'7 

(E2**Hz ) 

16xl0‘6 

7.7xlO'6 

2. OxlO"7 

26 

(E  -Hz) 

16xl0'6 

7.7xlO'6 

2. OxlO"7 

TABLE  2.3-2 

WHITE  NOISE  PARAMETERS 


BELL 

BASELINE 

BELL 

BALL-BEARING 

W1 

(E2/Hz) 

440 

650 

w2 

(E2/Hz ) 

440 

650 

w3 

(E2/Hz ) 

97 

300 

w4 

(E2/Hz ) 

97 

300 

w5 

(E2/Hz ) 

97 

1 

300 

?6 

(EZ/Hz) 

li 

i 

97 

300 

SIMULATION  RESULTS 


A  large  variety  of  simulation  results  were  obtained 
with  the  techniques  of  Chapter  2.  The  results  consist  of  grav 
ity  residuals  in  the  form  of  point  values  and  5  min,  15  mxn, 

1  deg  and  5  deg  means . 

In  order  to  examine  the  sensitivity  of  the  results  to 
the  local  characteristics  of  the  gravity  field,  three  differ¬ 
ent  field  models  were  used  in  the  simulations: 

•  Attentuated  White  Noise  (AWN)  gravity 

model 

•  Baseline  gravity  model 

•  Active  gravity  model. 

These  models  are  discussed  in  Section  3.1. 

The  simulation  results  are  given  in  Sections  3.2 

through  3.5.  Section  3.2  presents  results  for  nine  different 
survey  possibilities.  Sections  3.3,  3.4  and  3.5  examine  the 
sensitivity  of  specific  survey  alternatives  to  variations  in 
survey  parameters. 

3.1  GRAVITY  FIELD  MODELS 

A  stationary  gravity  field  model  is  completely  speci¬ 
fied  by  the  spectral  density  of  the  anomalous  surface  potent  it 
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o  ’  o 

through  the  frequency-domain  approach  given  in  Section  2.1. 
The  three  models  used  in  the  simulations  are  discussed  below. 

A  single  shell  of  the  Attenuated  White  Noise  model 
(Ref.  13)  contributes  a  surface  potential  with  a  spectral  den¬ 
sity  of  the  form 

4>(k)  (s)  =  87xd£ct£  e"47lDks  (3.1-1) 

This  spectral  density  can  be  viewed  as  arising  from  a  spherical 

shell  at  depth  below  the  surface  of  the  earth  on  which  the 

potential  is  white  (hence  the  name  Attenuated  White  Noise)  and 

2 

such  that  the  surface  potential  has  variance  cr£.  The  complete 
AWN  model  consists  of  five  independent  shells;  i.e., 

5 

4>t  ,(,)=£  *<k)  (3.1-2) 

°’  0  k=l 

Global  data  were  used  to  fit  the  ten  parameters  of  the  model 
in  Ref.  13.  The  resulting  values  of  the  parameters  are  given 
in  Table  3.1-1. 


TABLE  3.1-1 
AWN  MODEL  PARAMETERS 
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The  Baseline  and  Active  field  models  are  members  of 
the  same  class  of  models.  They  are  both  sums  of  two  third- 
order  Markov  models  (Ref.  14).  A  single  third-order  Markov 
model  has  a  surface  potential  with  a  spectral  density  of  the 
form 


4>(k)(s)  = 


10Mfrk 


<2ns>2]7/2 


(3.1-3) 


2 

where  1/0^  is  called  the  characteristic  distance  and  is  the 
variance  of  the  surface  potential.  For  the  Baseline  and  Active 
models 


2 

T  =  Y]  *(k)  (s)  (3.1-4) 

°*  °  k=l 

(k) 

with  4> v  7  as  in  Eq.  3.1-3.  The  four  parameters  of  the  Base¬ 
line  and  Active  models  were  obtained  by  fitting  to  data  in  the 
North  Atlantic  and  in  the  Bonin  Trench,  respectively.  Their 
values  are  given  in  Tables  3.1-2  and  3.1-3. 

TABLE  3.1-2 

BASELINE  MODEL  PARAMETERS 


K 

1/0  k 
(km) 

°k 

0  K  0 
(rn  /secz ) 

1 

27.78 

16.00 

2 

370.4 

91.43 

TABLE  3.1-3 

ACTIVE  MODEL  PARAMETERS 


K 

VPU 

(km) 

ak 

(mvsec^ ) 

1 

22.22 

29.98 

2 

350.03 

103.68 

Since  the  spectral  densities  given  in  Eqs.  3.1-1  and 
3.1-3  are  functions  of  the  frequency  magnitude  s,  the  three 
models  are  isotropic  .  Graphs  of  the  spectral  densities  of 
the  anomalous  surface  potential  for  the  AWN,  Baseline  and  Active 
models  as  functions  of  s  are  given  in  Fig.  3.1-1.  The  AWN  model 


1  1  1  ‘  1  ‘  4  L  — 1  I  “  1  I  .  1.  1  I  ■  I  ,  .  .  I  .  ,  !  J 

10.000  1000  100  10  1 

WAVELENGTH  |km) 

Figure  3.1-1  Spectral  Densities  of  the  Anomalous  Surface 

Potential  for  the  AWN,  Baseline,  and  Active 
Models 


*The  methodology  of  Chapter  2  is  applicable  to  non-isotropic 
models  as  well. 
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has  very  little  energy  at  high  frequencies;  it  represents  a  smooth 
slowly  varying  local  gravity  field.  On  the  other  hand,  the  Active 
model  contains  a  substantial  amount  of  high-frequency  energy  since 
as  its  name  indicates,  it  represents  a  very  active  local  gravity 
field.  The  Baseline  model  corresponds  to  an  intermediate  case. 

For  reference  purposes.  Table  3.1-4  presents  the  rms 
values  of  the  gravity  anomaly,  the  deflections  of  the  vertical 
and  gradients  of  the  anomalous  potential  at  the  surface  for 
the  three  models.  Shown  in  parentheses  are  the  1  deg  means. 


TABLE  3.1-4 

RMS  VALUES  OF  MODEL  QUANTITIES* 


QUANTITY 

AWN  MODEL 

BASELINE  MODEL 

ACTIVE  MODEL 

Gravity  Anomaly 

42.7 

51.1 

112.5 

(ragal) 

(36.2) 

(29.8) 

(48.1) 

East  Deflection  of  the  Vertical 

6.8 

7.6 

16.8 

(sec) 

(5.4) 

(4.4) 

(7.2) 

North  Deflection  of  the  Vertical 

6.8 

7.6 

16.8 

(sec) 

(5.4) 

(4.4) 

(7.2) 

East-East  Gradient 

13.8 

20.7 

60.6 

(E) 

(3.5) 

(3.8) 

(8.1) 

North-North  Gradient 

13.8 

20.7 

60.6 

(E) 

(3.5) 

(3.8) 

(8.1) 

Vertical-Vertical  Gradient 

22.6 

33.8 

99.0 

(E) 

(5.7) 

(6.2) 

(13.0) 

East-North  Gradient 

8.0 

12.0 

35.0 

(E) 

(2.0) 

(2.2) 

(4.4) 

East-Vertical  Gradient 

16.0 

23.9 

70.0 

(E) 

(4.0) 

(4.4) 

(9.2) 

North-Vertical  Gradient 

16.0 

23.9 

70.0 

(E) 

(4.0) 

(4.4) 

(9.2) 

*1  deg  means  shown  in  parentheses. 
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COMPARISON  OF  SURVEY  ALTERNATIVES 


To  illustrate  the  versatility  of  the  methodology  de¬ 
scribed  in  Chapter  2,  nine  survey  possibilities  were  simulated 
using  the  three  field  models  introduced  in  the  previous  section. 
Post-survey  rms  residuals  in  the  deflections  of  the  vertical  and 
the  gravity  anomaly  are  presented  in  Tables  3.2-1,  3.2-2,  and 
3.2-3  for  the  AWN,  Baseline,  and  Active  models,  respectively. 
All  results  correspond  to  a  survey  in  an  equatorial  region. 

Each  table  contains  10  columns.  The  first  column, 
labeled  NONE,  corresponds  to  the  unsurveyed  field.  The  acro¬ 
nyms  ALT,  SST,  SHIP,  GRAV  and  GRAD,  used  in  the  remaining  nine 
columns,  represent  satellite  radar  altimetry,  satellite- to- 
satellite  tracking,  ship  gravimetry,  land-based  gravimetry  and 
airborne  gradiometry,  respectively.  Each  column  corresponds 
to  the  survey  combination  indicated  by  the  acronyms  in  its 
header.  A  description  of  the  survey  parameters  is  given  next. 

Satellite  Radar-Altimeter  Survey  -  The  GEOS -3  geometry 
and  parameters  were  used.  The  equatorial  separation  of  GEOS-3 
groundtracks  was  taken  as  30  nm. 

Satellite- to-Satellite  Tracking  -  The  low  satellite's 
altitude  is  150  km.  Its  orbit  inclination  is  94  deg.  Measure¬ 
ments  of  line-of-sight  range-rate  are  taken  every  10  sec  and 
the  noise  level  has  a  standard  deviation  of  1  pm/sec  per  meas¬ 
urement.  The  separation  of  the  groundtracks  of  the  low  satel¬ 
lite  is  30  nm  at  the  equator.  The  high  satellite  is  directly 
above  the  region  where  estimates  are  sought. 

Shipborne  Gravimetry  -  Spacings  between  contiguous 
survey  tracks  were  selected  following  NAVOCEANO's  methodology. 
The  remaining  parameters  are  defined  in  Subsection  2.3.3. 
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TABLE  3.2-1 

AWN  MODEL  SIMULATION  RESULTS 


SURVEY  TYPE 

NONE 

SST 

ALT 

SST 

ALT 

GRAV 

GRAD 

ALT 

SHIP 

SST 

GRAD 

ALT 

GRAD 

SST 

ALT 

GRAD 

H 

RMS  North  Deflection 
of  the  Vertical*  (sec) 

6.8 

(5.4) 

2.1 

(0.5) 

1.4 

(0.4) 

1.3 

(0.1) 

0.8 

0.7 

0.5 

a 

0.4 

RMS  East  Deflection 
of  the  Vertical*  (sec) 

6.8 

(5.4) 

1.9 

(0.5) 

1.5 

(0.5) 

1.4 

(0.1) 

0.6 

0.4 

0.4 

0.4 

0.3 

RMS  Gravity  Anomaly* 
(mgal) 

42.7 

(36.2) 

13.2 

(3.2) 

9.6 

(2.8) 

9.1 

(0.7) 

■ 

a 

3.0 

3.0 

2.8 

2.4 

*1  deg  means  shown  in  parenthesis. 


TABLE  3.2-2 

BASELINE  MODEL  SIMULATION  RESULTS 


SURVEY  TYPE 


NONE 

SST 

ALT 

SST 

ALT 

GRAV 

GRAD 

ALT 

SHIP 

SST 

GRAD 

ALT 

GRAD 

SST 

ALT 

GRAD 

ALT 

GRAD 

SHIP 

RMS  North  Deflection 
of  the  Vertical*  (sec) 

7.6 

(4.4) 

4.5 

(0.5) 

M 

2.1 

(0.2) 

1.0 

0.6 

0.6 

0.5 

n 

0.4 

RMS  East  Deflection 
of  the  Vertical*  (5e?) 

7.6 

(4.4) 

4.2 

(0.4) 

2.8 

(0.5) 

2.8 

(0.3) 

0.6 

0.6 

0.4 

0.4 

0.4 

0.3 

RMS  Gravity  Anomaly* 
(mgal) 

51.1 

(29.8) 

29.2 

(3.0) 

16.9 

(3.2) 

16.7 

(1.7) 

5.2 

4.2 

3.4 

3.2 

3.0 

2.5 

*1  deg  means  shown  in  parenthesis. 
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TABLE  3.2-3 

ACTIVE  MODEL  SIMULATION  RESULTS 


SURVEY  TYPE 

NONE 

SST 

ALT 

SST 

GRAV 

ALT 

SST 

ALT 

SST 

PH 

ALT 

GRAD 

SHIP 

GRAD 

GRAD 

ALT 

GRAD 

m 

RMS  North  Deflection 

16.8 

12.5 

5.4 

5.4 

1.7 

■ 

1.2 

of  the  Vertical*  (sec) 

(7.2) 

(1.6) 

(0.6) 

(0.4) 

O .  y 

0 .  y 

U .  o 

RMS  East  Deflection 

16.8 

11.4 

8.0 

7.9 

1.1 

i 

0.8 

0.8 

0.8 

of  the  Vertical*  (sec) 

(7.2) 

(1.6) 

(1.0) 

(0.8) 

B 

RMS  Gravity  Anomaly* 

112.5 

80.2 

46.0 

45.7 

■ 

(■gal ) 

(48.1) 

(10.5) 

(5.5) 

(4.4) 

m 

m 

*1  deg  means  shown  in  parenthesis. 


Land-Based  Gravimetry  -  All  parameters  are  as  described 
in  Subsection  2.3.3. 

Airborne  Gradiometry  -  Data  are  assumed  to  be  collected 
with  the  Bell  Baseline  gradiometer  triad  aboard  an  aircraft  fly¬ 
ing  at  a  speed  of  300  kt  at  an  altitude  of  20,000  ft.  Measure¬ 
ments  are  taken  every  10  seconds  on  parallel  tracks  spaced  10  km 
apart . 

Three  of  the  five  sensors  considered  provide  informa¬ 
tion  on  the  long-wavelength  content  of  the  gravity  field. 

They  are  SST,  ALT,  and  GRAV.  The  other  two  sensors,  GRAD  and 
SHIP,  recover  the  short  wavelengths  of  the  gravity  field.  In 
the  case  of  the  gravimetric  surveys  (GRAV  and  SHIP),  this  dif¬ 
ference  is  primarily  a  function  of  the  density  of  coverage. 

For  the  other  sensors  (SST,  ALT,  and  GRAD),  the  bandwidth  of 
recovery  is  inherent  in  the  physical  characteristics  of  the 


sensor . 


In  Tables  3.2-1,  3.2-2,  and  3.2-3  all  columns  contain¬ 
ing  ALT  or  SHIP  in  the  header  correspond  to  ocean  surveys. 

The  column  headed  by  GRAV/GRAD  represents  a  land  survey.  The 
two  columns  SST  and  SST/GRAD  apply  to  ocean  and  to  overland 
surveys . 

The  three  tables  present  the  point  value  rms  residuals 
in  the  north  deflection  of  the  vertical,  the  east  deflection 
of  the  vertical,  and  the  gravity  anomaly.  The  numbers  in  paren 
thesis  are  the  rms  values  of  the  residual  1  deg  means.  These 
are  given  only  for  those  surveys  in  which  long-wavelength  sen¬ 
sors  are  used  exclusively. 

The  results  given  in  Tables  3.2-1,  3.2-2,  and  3.2-3 
are  discussed  in  detail  below.  First,  consider  the  order  of 
the  columns  in  the  tables.  Note  that  the  various  columns  in 
the  three  tables  are  ordered  in  the  same  sequence.  This  order 
corresponds  to  decreasing  rms  values  for  the  residual  point 
gravity  anomaly  in  the  Baseline  and  AWN  models.  In  the  case 
of  the  Active  gravity  model,  the  gravity  anomaly  residual  for 
the  combination  ALT/SHIP  appears  out  of  sequence.  The  reason 
is  that  the  rms  residuals  shown  in  the  tables  correspond  to 
the  actual  ship  survey  that  NAVOCEANO  would  conduct.  The  cri¬ 
terion  used  by  NAVOCEANO  to  select  the  spacing  between  con¬ 
tiguous  ship  tracks  would  choose  the  same  spacing  for  the  AWN 
and  Baseline  models,  but  it  would  choose  a  denser  collection 
of  tracks  for  a  region  described  by  the  Active  gravity  model. 

If  the  same  ship  track  spacing  is  used  for  all  models,  the 
results  for  the  Active  model  appear  in  the  same  order  as  those 
of  the  other  two  models. 

For  most  surveys,  the  rms  values  of  the  residuals  in 
the  north  and  in  the  east  deflections  of  the  vertical  are  not 
the  same.  This  illustrates  the  fact  that  the  post-survey  gravi 
ty  residuals  are  not  isotropic.  Even  though  the  three  models 
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used  in  the  simulations  are  isotropic,  the  gravity  residuals 
turn  out  to  be  anisotropic  for  all  surveys.  The  primary  reason 
for  the  anisotropy  is  the  geometry  of  the  surveys.  For  example, 
for  an  airborne  gradiometer  survey  the  survey  tracks  are  east- 
west  and  samples  along  a  track  are  spaced  approximately  1.5  km 
apart,  while  in  the  cross-track  direction  (north-south)  the 
distance  between  samples  is  10  km.  In  the  case  of  the  SST 
survey,  the  rms  of  the  residual  east  deflection  of  the  verti¬ 
cal  is  smaller  than  that  of  the  north  deflection  of  the  verti¬ 
cal  because  the  10  sec  range-rate  sampling  implies  a  larger 
spacing  between  consecutive  samples  on  a  groundtrack  (roughly 
north-south)  than  the  spacing  between  continuous  groundtracks 
(east-west) . 

The  advantages  of  multisensor  surveys  are  clear  from 
the  tables.  The  smallest  residuals  are  obtained  when  a  combi¬ 
nation  of  sensors  is  used  to  cover  the  entire  range  of  fre¬ 
quencies.  This  involves  at  least  one  long-wavelength  sensor 
and  one  short-wavelength  sensor.  Little  is  gained  by  combin¬ 
ing  two  long-wavelength  sensors,  as  can  be  seen  by  comparing 
the  columns  SST  and  ALT  with  the  column  headed  by  SST/ALT  in 
each  of  the  tables.  The  best  two-sensor  results  are  obtained 
with  the  ALT/GRAD  sensor  combination  in  the  cases  of  the  AWN 
and  the  Baseline  gravity  models.  For  the  Active  model,  the 
ALT/SHIP  sensor  combination  provides  the  best  recovery  using 
only  two  sensors.  In  each  case,  the  best  two-sensor  result  is 
not  very  much  improved  when  a  third  sensor  is  added. 

3.3  SATELLITE  ALTIMETRY  -  SENSITIVITY  TO  TRACK  SPACING 

This  section  presents  a  study  of  the  sensitivity  of  the 
gravity  recovery  to  satellite  altimeter  track  spacing.  Both 
GEOS-3  and  SEASAT-1  altimeters  are  considered.  Some  important 
implications  of  the  survey  geometry  are  discerned  and  discussed. 


Figures  3.3-1,  3.3-2,  and  3.3-3,  respectively,  show 
the  variation  with  equatorial  track  spacing  of  the  residual 
rms  gravity  anomaly,  east  deflection  of  the  vertical,  and  north 
deflection  of  the  vertical.  Results  for  the  AWN  and  Active 
models  for  recovery  at  the  equator  are  given  in  these  figures. 

The  horizontal  axes  correspond  to  the  equatorial  separation 
between  contiguous  tracks. 

At  first,  the  results  appear  surprising.  The  SEASAT-1 
altimeter  is  of  better  quality  than  the  GEOS-3  altimeter.  In 
addition,  SEASAT-1  has  a  larger  orbit  inclination  than  GEOS-3. 
Consequently,  it  is  expected  that  the  gravity  anomaly  and  the 
north  deflection  of  the  vertical  are  better  recovered  by  SEASAT-1 
than  by  GEOS-3  at  all  track  spacings  .  However,  Figs.  3.3-1 
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Figure  3.3-1  Sensitivity  of  Residual  Gravity  Anomaly 

to  Altimeter  Survey  Track  Spacing 


*For  the  east  deflection  of  the  vertical  (Fig.  3.3-2)  the  results 
are  expected  because  of  the  difference  in  orbit  inclination  be¬ 
tween  the  two  satellites. 
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Figure  3.3-3 
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EQUATORIAL  TRACK  SPACING  (nm) 

Sensitivity  of  Residual  North  Deflection  of  the 
Vertical  to  Altimeter  Survey  Track  Spacing 
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and  3.3-3  indicate  that  for  a  wide  range  of  track  spacings  the 
residuals  in  the  gravity  anomaly  and  the  north  deflection  of 
the  vertical  are  larger  for  SEASAT-1  than  for  GEOS-3.  Only 
for  very  short  track  spacings  do  SEASAT-1  results  show  an  im¬ 
provement  over  GEOS-3. 

The  origin  of  this  surprising  behavior  is  in  the  geome 
try  of  the  survey.  Figure  3.3-4  illustrates  the  geometric 
differences  between  SEASAT-1  and  GEOS-3  surveys.  For  the  same 
equatorial  separation,  any  two  parallel  SEASAT-1  tracks  are 
separated  by  a  larger  distance  in  the  north  direction  than  two 
GEOS-3  tracks.  More  fundamentally,  as  Fig.  3.3-5  illustrates, 
for  the  same  equatorial  separation  there  are  more  data  per 
unit  area  in  a  GEOS-3  survey  than  in  a  SEASAT-1  survey.  Thus, 
the  results  reflect  a  trade-off  between  data  quality  and  data 
quantity. 


R — 47720 


Figure  3.3-4  Geometric  Differences  Between  SEASAT-1 

and  GEOS-3  Surveys 
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Figure  3.3-5  Differences  in  Density  of  Coverage 

Between  SEASAT-1  and  GEOS -3  Surveys 


When  a  SEASAT-1  quality  altimeter  is  used  in  a  satel¬ 
lite  at  the  same  orbit  of  GE0S-3,  the  results  of  Figs.  3.3-6, 
3.3-7  and  3.3-8  are  obtained.  For  ease  in  referencing,  this 
case  has  been  labeled  SEASAT-l-A.  It  can  be  seen  from 
Figs.  3.3-6,  3.3-7  and  3.3-8  that,  as  expected,  SEASAT-l-A 
yields  consistently  better  recovery  than  GE0S-3  but  the  im¬ 
provement  is  evident  only  for  track  spacings  shorter  than 
40  nm. 

Residual  5  nun,  15  min,  1  deg,  and  5  deg  means  of  the 
gravity  anomaly  are  presented  in  Figs.  3.3-9  through  3.3-12. 
The  effects  of  the  geometric  differences  between  the  surveys 
are  clearly  seen  in  these  figures.  The  effects  of  the  differ¬ 
ences  between  the  quality  of  the  altimeters  corresponding  to 
the  GEOS-3  and  SEASAT  surveys  are  only  noticeable  in  the  re¬ 
covery  of  the  5  and  15  nun  means.  The  values  of  the  5  nun, 
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Figure  3.3-6  Sensitivity  of  Residual  Gravity  Anomaly  to 

Altimeter  Survey  Track  Spacing 
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Figure  3.3-7  Sensitivity  of  Residual  East  Deflection  of  the 

Vertical  to  Altimeter  Survey  Track  Spacing 
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3.3-8  Sensitivity  of  Residual  North  Deflection  of  the 
Vertical  to  Altimeter  Survey  Track  Spacing 


Figure  3.3-9  Satellite  Altimeter  Survey  Residual 

5  mTn  Means 
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Figure  3.3-10  Satellite  Altimeter  Survey  Residual 

15  min  Means 


Figure  3.3-11  Satellite  Altimeter  Survey  Residual 

1  deg  Means 
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EQUATORIAL  TRACK  SPACING  (nm) 

Figure  3.3-12  Satellite  Altimeter  Survey  Residual 

5  deg  Means 


15  min,  1  deg,  and  5  deg  means  for  the  unsurveyed  AWN  and  Active 
fields  are  given  as  a  reference  in  Table  3.3-1.  Note  that  the 
5  deg  means  of  the  unsurveyed  field  are  larger  for  the  AWN 
model  than  for  the  Active  model  because  the  AWN  model  has  more 
power  in  the  long  wavelength  portion  of  the  spectrum  than  the 
Active  model  (see  Fig.  3.1-1).  This  explains  why  the  graphs 
of  Fig.  3.3-12  cross  for  sparse  surveys. 

TABLE  3.3-1 

SPATIALLY  AVERAGED  UNSURVEYED  GRAVITY  ANOMALY 


5  min 
(mgal) 


15  rain 
(mgal) 


1  deg 
(mgal) 


41.6 


40.2 


36.2 


26.1 


Active 


110.0 


97.7 


48.1 


20.1 


Figures  3.3-11  and  3.3-12  indicate  that  there  are 
fundamental  lower  bounds  on  the  recovery  of  the  1  deg  and  5  deg 
means  of  the  gravity  anomaly  from  satellite  altimetry  data. 
These  bounds  (approximately  2.6  mgal  for  the  1  deg  means  and 
1.2  mgal  for  the  5  deg  means)  are  independent  of  the  gravity 
field  model  used  in  the  analysis  and  are  also  independent  of 
the  survey  geometry  and  the  quality  of  the  altimeter  used. 

3.4  SST  -  SENSITIVITY  TO  LOW  SATELLITE'S  ALTITUDE 

A  study  was  conducted  of  the  sensitivity  of  SST  sur¬ 
vey  residual  gravity  errors  with  respect  to  the  height  of  the 
low  satellite.  Measurements  of  line-of-sight  range-rate  were 
considered  to  be  taken  every  10  sec  and  to  have  an  uncertainty 
with  a  standard  deviation  of  1  pm/sec .  The  inclination  of  the 
orbit  of  the  low  satellite  was  taken  as  94  deg  and  the  equa¬ 
torial  separation  between  neighboring  groundtracks  was  kept 
constant  at  30  run.  The  altitude  of  the  low  satellite  over  the 
earth's  surface  was  allowed  to  vary  between  100  km  and  800  km. 

Results  for  gravity  recovery  from  an  equatorial  region 
directly  under  the  high  satellite  were  obtained  for  the  AWN 
and  Active  models.  The  point  residuals  in  the  gravity  anomaly 
are  presented  in  Fig.  3.4-1.  Figures  3.4-2  through  3.4-5  show 
the  5  min,  15  min,  1  deg,  and  5  deg  means  of  the  residuals  in 
the  gravity  anomaly. 

As  the  height  of  the  low  satellite  decreases,  the 
spacing  between  consecutive  samples  along  a  groundtrack  in- 
creases  according  to  Eq.  2.3-1  .  On  the  other  hand,  the  closer 

*Along- track  sample  spacing  for  altitudes  of  100  km  and  800  km 
are  66.2  km  and  77.2  km,  respectively. 
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ALTITUDE  OF  LOW  SATELLITE  (km) 

Figure  3.4-1  SST  Survey  RMS  Residual  Point 

Gravity  Anomaly 
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Figure  3.4-2  SST^Survey  Gravity  Anomaly  Residual 

5  min  Means 
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Figure  3.4-3  SST^Survey  Gravity  Anomaly  Residual 

15  min  Means 
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Figure  3.4-4  SST  Survey  Gravity  Anomaly  Residual 

1  deg  Means 
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Figure  3.4-5  SST  Survey  Gravity  Anomaly  Residual 

5  deg  Means 

the  low  satellite  is  to  the  earth's  surface,  the  stronger  the 
anomalous  gravity  field  and  the  better  the  signal- to-noise 
ratio  of  the  measurements.  By  far,  the  dominant  effect  is  the 
quality  of  the  data.  Thus,  even  though  the  data  become  sparser 
as  the  height  of  the  low  satellite  decreases,  residual  gravity 
errors  are  reduced. 

Table  3.4-1  gives  the  residual  5  min,  15  min,  1  deg, 
and  5  deg  means  of  the  gravity  anomaly  as  fractions  of  the 
corresponding  unsurveyed  field  values  for  low  satellite's  alti¬ 
tudes  of  100  km  and  800  km.  The  fractions  have  been  rounded- 
off  to  the  nearest  1/100.  The  5  deg  means  are  very  well  re¬ 
covered  from  SST  data.  In  general,  the  fractional  recovery  of 
the  gravity  anomaly  is  better  for  the  AWN  model  than  for  the 
Active  model  because  most  of  the  energy  in  the  field  repre¬ 
sented  by  the  AWN  model  is  contained  in  the  long-wavelength 
portion  of  the  spectrum  while  in  the  case  of  the  Active  model 
there  is  a  large  amount  of  energy  at  high  frequencies  which  is 
unobservable  in  the  data. 


TABLE  3.4-1 

RESIDUAL  GRAVITY  ANOMALY  AS  A  FRACTION  OF 
THE  UNSURVEYED  FIELD  VALUES 


HEIGHT  OF 

LOW  SATELLITE 

MODEL 

5  min 
MEAN 

15  min 
MEAN 

1  deg 
MEAN 

5  deg 
MEAN 

100  km 

AWN 

0.28 

0.22 

0.04 

0.01 

100  km 

Active 

0.68 

0.60 

0.17 

0.05 

800  km 

AWN 

0.85 

0.84 

0.64 

0.09 

800  km 

Active 

0.97 

0.96 

0.83 

0.17 

3.5  AIRBORNE  GRADIOMETRY  -  SENSITIVITY  TO  TRACK  SPACING 

This  section  presents  results  on  the  sensitivity  of 
residual  gravity  errors  from  a  combined  SST/Airborne  gradio- 
metric  survey.  The  gradiometer  instruments  considered  are 
the  Bell  Ball-Bearing  and  the  CSDL  triads  discussed  in 
Subsection  2.3.4. 

For  the  SST  survey,  the  low  satellite's  altitude 
above  the  earth's  surface  was  taken  as  120  km.  The  rest  of 
the  parameters  were  the  same  as  those  used  in  obtaining  the 
results  of  the  previous  section.  The  gradiometer  survey  was 
assumed  to  be  conducted  at  an  altitude  of  20,000  ft  and  at  a 
speed  of  300  kt.  The  interval  between  successive  gradiometer 
samples  was  taken  as  10  sec.  Track  spacing  was  allowed  to 
vary  between  10  km  and  100  km. 

The  results  for  the  rms  of  the  residual  gravity  anom¬ 
aly  are  presented  in  Fig.  3.5-1.  Note  that  the  differences  in 
gravity  recovery  between  the  Bell  and  CSDL  instruments  are 
exaggerated  in  Fig.  3.5-1  because,  as  noted  in  Subsection  2.3.4, 
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Figure  3.5-1  Sensitivity  of  SST/Airborne  Gradiometer  Survey 

Residual  Gravity  Anomaly  to  Gradiometer  Survey 
Track  Spacing 

the  error  model  for  the  CSDL  instrument  does  not  include  the 
effects  of  translational  and  rotational  vibration* ,  while  the 
Bell  Ball-Bearing  instrument  error  model  does  include  these 
effects . 


There  are  three  possibilities  to  further  decrease  the 
residual  gravity  anomaly.  First ,  the  slope  of  the  residual 
gravity  anomaly  curve  at  the  lower  end  of  the  track  spacing 
scale  in  Fig.  3.5-1  indicates  that  the  full  benefit  of  an  air¬ 
borne  gradiometer  survey  has  not  been  attained  at  a  track  spac 
ing  of  10  km.  Second .  residual  gravity  anomaly  can  be  reduced 
by  lowering  the  height  of  the  gradiometer  survey  (20,000  ft). 
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Third,  the  east-west  gradiometer  survey  could  be  complemented 
with  a  north-south  gradiometer  survey  to  provide  short- 
wavelength  recovery  over  the  frequency-domain  region  where  the 
red  noise  in  the  east-west  gradiometer  survey  significantly 
corrupts  the  measurements. 
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4. 


CONCLUSIONS  AND  RECOMMENDATIONS 


4.1  CONCLUSIONS 

A  new  method  and  a  computer  program  have  been  devel¬ 
oped  for  the  analysis  of  multisensor  surveys  of  the  gravity 
field.  The  method  is  ideally  suited  for  the  design  of  multi¬ 
sensor  surveys  to  achieve  a  desired  level  of  gravity  recovery 
accuracy. 

The  approach,  which  utilizes  a  flat-earth  approxima¬ 
tion  in  the  survey  region,  is  based  on  an  extension  of  optimal 
Wiener  smoothing  theory.  The  Fourier  transformation  permits 
the  efficient  evaluation  of  the  statistics  of  the  post-survey 
residuals  in  the  form  of  their  average  spectral  density.  Root- 
mean-square  (rms)  values,  covariances  and  crosscovariances  of 
point  and  spatial  averages  of  the  residuals  are  readily  com¬ 
puted  from  their  average  spectral  density  through  numerical 
integration. 

Survey  error  models  and  gravity  recovery  simulation 
results  were  given  for  a  variety  of  combinations  of 

•  Satellite  radar  altimeters 

•  Satellite-to-satellite  tracking  in  a 
high-low  configuration 

•  Land-based/shipborne  gravimetry 

•  Airborne  gradioraetry. 

The  examples  given  illustrate  the  versatility  of  the  methodology. 
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RECOMMENDATIONS 


There  are  several  important  areas  of  investigation 
for  future  extensions  of  the  methodology: 


•  Effects  of  survey  irregularities  -  The 
effects  of  gaps  in  satellite  data  and 
the  seemingly  random  measurement  pattern 
associated  with  land-based  gravimetric 
surveys  would  be  incorporated  in  the 
analysis . 

•  Additional  sensors  -  Other  survey  possi¬ 
bilities  such  as  satellite- to-satellite 
tracking  in  a  low-low  configuration  or 
inertial  gravimetry  can  be  included  in 
the  formulation. 


APPENDIX  A 

THE  FOURIER  TRANSFORMATION 


This  appendix  presents  a  compilation  of  those  defini¬ 
tions  and  results  of  Fourier  analysis  employed  in  this  report. 
This  appendix  is  not  intended  as  a  treatise  on  the  subject  of 
Fourier  analysis.  The  material  is  presented  as  an  aid  in  under¬ 
standing  the  application  of  the  concepts  rather  than  as  an 
investigation  of  the  mathematical  properties  of  the  concepts 
themselves.  The  reader  interested  in  a  more  complete  presen¬ 
tation  is  referred  to  various  excellent  books  on  the  subject 
such  as  Refs.  3,  15,  16,  and  17. 


A. 1  DEFINITIONS  OF  FOURIER  TRANSFORMS 

Let  g(x)  be  a  real-valued  function  defined  for  all 
T 

points  x  =  (x^,X2>  of  a  cartesian  coordinate  system.  The 
Fourier  transform  of  g,  g,  is  a  complex-valued  function  defined 
by 


00 

g(s)  =  JJ  g(x)e“l2n<-,->  dx-^dx. 


( A- 1 ) 


T 

where  i  =  /-I ,  s  =  (s^,S2)  >  and  where  <x,s>  is  the  scalar 
product  of  the  vectors  x  and  s  given  by 


T 

<X,S>  =  X  S  =  X1S1  +  X2S2 


(A-2 ) 


The  vector  s  is  referred  to  as  planar  frequency.  Its  com¬ 
ponents,  s^  and  S2»  are  real  numbers  which,  as  shown  below, 
correspond  to  physical  frequencies  measured  in  the  directions 


A- 1 


of  the  axes  of  the  coordinate  system.  The  fact  that  g  is 
complex-valued  is  solely  due  to  the  appearance  of  the  factor  i 
in  the  exponential  of  the  integrand  in  Eq.  A-l.  Note  that 
since  g(x)  is  real-valued,^ 


g(s)  =  (g(-s)] 


(A-3) 


The  inverse  Fourier  transformation  permits  the  re¬ 
covery  of  g(x)  from  g(s): 


CD 

g(x)  =  JJ g(s)el2n<-’->  ds1ds2 


(A-4) 


The  formulas  for  the  Fourier  transform  and  the  inverse  Fourier 
transform  are  completely  symmetric  except  for  a  change  of  sign 
in  the  exponential.  The  functions  g  and  g  are  called  a  Fourier 
transform  pair.  Given  any  one  of  them,  the  other  is  uniquely 
determined  within  an  equivalence  class  (Ref.  16). 


In  order  to  give  a  physical  interpretation  of  the 
Fourier  transformation,  it  is  convenient  to  decompose  g(s) 
into  its  real  and  imaginary  parts.  Let 


g(s)  =  gr(s)  +  ig^ (s)  (A-5 ) 

where  gr  and  g^  are  real.  From  Eq.  A-3  it  follows  that 

gr(s)  -  gr(-s)  (A-6) 

and 

gf(s)  =  -gi(-s)  (A-7 ) 

tA  superscript  asterisk  denotes  complex  conjugation  when  attached 
to  a  scalar  and  conjugate  transpose  when  attached  to  a  matrix  or 
a  vector. 
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Using  the  identity 

i2n<x,s>  0  „ 

e  -  =  cos  2n<x,s>  +  1  sxn  2n<x,s>  (A-8) 

and  Eq.  A-5,  the  inverse  Fourier  transformation  can  be  written 
as 


=  J'J'  [fL r(s)  cos27i<x,s>  -  g^(s)  sin2n<x , s>  Jds^ds. 
-00 


00 

+  i  fj  lgr(s)  sin2n<x,s>  +  gi(s)  cos2n<x , s>  ]ds1ds2 
-00 


( A-9 ) 


From  Eqs.  A-2,  A-6  and  A-7  it  follows  that  the  second  integral 
in  the  right-hand  side  of  Eq.  A-9  vanishes  and  that 

00  00 

■*>■//  All<sl  ,S2>  cos2ns^x^  cos2ns2X2  ds^ds2 


o  o 


00  00 


If  A22(s  ^,S2)  sin2;tS2Xj,  sin2ns2X2  ds^ds2 


o  o 


00  00 

//  A12(s1,s2)  cos2ns2X^  sin2ns2X2  ds^ds. 


o  o 


00  00 

J *  J"  A22(S2»S2)  sin2ns2X2  cos2ns2X2  ds2ds2 


o  o 


A-  3 


( A- 10 ) 


where 


A11 ( S1 ’ s2  >  =  2lSr(s1,s2)  +  gr(s1,-s2)l 


A22(s1,s2)  =  2(-gr(slfs2)  +  gr(s1,-s2)] 


Aj 2 ( s i , s 2 )  2  [  -g^  ( s i  , s0  )  + 


( A-ll ) 


®iv  1*  2' 


A21(s1,s2)  =  2[-gi(s1,s2)  -  gi(s1,-s2)] 


Equation  A-10  is  a  representation  of  g(x)  as  a  sum 
over  frequency  of  products  of  trigonometric  functions  for  the 
four  possible  phase  combinations.  Only  positive  frequencies 
appear  in  Eq.  A-10.  The  Fourier  transform  g(s)  associates 
with  each  pair  of  positive  frequencies  (s^,s2)  eight  real  num¬ 
bers  corresponding  to  the  real  and  imaginary  parts  of  g(s)  fo 
the  four  vector  planar  frequencies  s  =  (s^,s2>  ,  s  =  (s^,-s2) 
s  =  (-s^,s2)  and  s  =  (-s^,-s2)^.  Of  these  eight  numbers  only 
four  are  independent  (Eq.  A-3)  and  the  relation  between  these 
numbers  and  the  coefficients  of  the  expansion  A-10  is  given  by 
A-ll.  The  definition  of  the  Fourier  transform  (Eq.  A-l)  as  a 
complex-valued  function  of  positive  and  negative  frequencies 
is  a  compact  mathematically  convenient  way  of  representing  the 
four  amplitudes  of  the  trigonometric  expansion  A-10.  The  in¬ 
version  formula  (Eq.  A-4)  is  entirely  equivalent  to  the  expan¬ 
sion  A-10. 


Classically,  the  most  important  application  of  the 
Fourier  transformation  has  been  to  the  solution  of  differen¬ 
tial  equations.  This  application  is  a  consequence  of  the  fol¬ 
lowing  property  of  Fourier  transforms:  Consider  the  n-th 
partial  derivative  of  g(x)  with  respect  to  x^  (k=l,2).  From 
the  formula  for  the  inverse  Fourier  transform  (Eq.  A-4) 
3ng(x)/3xP  can  be  written  as 


H  '1 
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bl  *<5>  =JJ  -n  l8(s)  dSlds2  (A-12) 

3xk  -£  3xk 

With  the  help  of  Eq.  A-2,  the  partial  derivative  of  the  inte¬ 
grand  can  be  evaluated  to  yield 


.n 


9x 


n 


g(x) 


00 

■// 


[<i2nsk)n  g(s)] 


ei<x,s>  d  ^  (A-13) 


from  which  it  follows  that  the  Fourier  transform  of  ang(x)/3x 
is  the  Fourier  transform  of  g  multiplied  by  (i2nsk)n. 


As  an  example  of  the  application  of  this  property  of 
the  Fourier  transformation,  consider  the  solution  of  Laplace's 
equation  with  a  boundary  condition  on  a  plane.  Let  T  (x)  be 

T 

the  potential  at  height  z  above  the  point  x  =  (x^,X2>  .  At 
z=0,  the  potential  is  known  and  is  given  by  Tq(x).  For  z>0, 

T  (x)  satisfies 

Z 

2  2  2 

T  <x)  ♦  T  (x)  +  T  (x)  =  0  <A-14) 

axz  z  ax2  azz  z 

In  addition,  it  is  required  that  T  (x)  remain  bounded  as  z  -*■  ® 

Z  ” " 

A 

Let  T  (s)  be  the  Fourier  transform  of  the  potential 

Z  ““ 

on  the  plane  at  height  z.  Applying  the  Fourier  transformation 
to  Eq.  A-14  results  in 


-4nzsfT  (s) 


4n2s?T  ( s) 


=  0 


(A-15 ) 


c  JC 


from  which  it  follows  that 


A  A  A  A 

Tz(s)  =  4nVTz<s) 


(A-16) 


where 

s  =  <s2  +  s^)1^2  (A-17) 

Equation  A-16  is  an  ordinary  differential  equation  for  each 
value  of  s.  Its  solution  is  of  the  form 

T  (s)  =  A  e“2nsz  +  B  e+2nsz  (A-18) 

where  A  and  B  are  to  be  determined.  From  the  condition  .that 
the  potential  must  be  bounded  as  z  +  ®  it  follows  that  B=0. 

A 

The  value  of  A  is  obtained  by  setting  z=0  to  yield  A  =  Tq(s). 
Thus,  the  solution  to  Eq.  A-16  is 

Tz(s)  =  e'2nsz  Tq(s)  (A-19) 

This  equation  states  that  the  Fourier  transforms  of  the  poten¬ 
tial  at  height  z  and  at  height  0  are  simply  related  by  the 

r\ 

attenuation  factor  e  nsz .  In  the  space  domain,  the  solution 

is  obtained  by  using  the  inverse  Fourier  transform;  i.e., 

00 

Tz(x)  =  JJ  TQ(s)  e‘2nsz  ei2n<x,s>  dSi<is2  (A-20) 

-00 

In  the  discrete  case  the  Fourier  transform  is  called 
the  finite  Fourier  transform  and  is  defined  in  a  manner  analo¬ 
gous  to  Eq.  A-l.  Let  M  be  an  orthogonal  grid  of  points  on  the 
plane  as  in  Fig.  A-l.  The  normalized  finite  Fourier  transform 


with  period  1  in  each  direction;  i.e.,  if  /^(A.m)1  is  a  vector 
with  arbitrary  integer-valued  entries  then 


g<£+A)  =  £  G<°>  e‘i2"<a’E+-> 
Q 

=  G(£) 


A-22) 


because 

e-i2n<Q  ,j>+A>  _  e-i2n<ft,j>>  .  e-i2n<Q,A> 


(A-23) 


and 


t-i2n<fl,A>  _  e-i2n( jA+km)  _  ^ 


(A-24) 


for  all  possible  values  of  j ,  A,  k  and  m. 

The  inverse  of  the  normalized  finite  Fourier  trans¬ 
form  is  given  by 

h  %  . 

G(Q)  =  f  f  G(fi)  ei2/l<-’^>  dPl  dp2  (A-25) 

~h  ~h 

The  limits  of  integration  in  Eq.  A-25  can  be  modified  as  long 
as  they  span  a  full  period  of  G(j>)  in  each  dimension. 

The  vector  p  =  (p^,P2)  in  ^qs.  A-21  and  A-25  is  called 
normalized  planar  frequency  because  it  does  not  carry  any  infor¬ 
mation  as  to  the  physical  wavelengths  implied  by  p^  and  p2< 

The  reason  is  that  the  function  G  has  been  viewed  as  being 
defined  on  the  integer  coordinates  of  the  grid  and  the  spacing 
between  consecutive  points  in  each  direction  is  one  unit  ir¬ 
respective  of  the  physical  distance  involved. 


(j ,k)T  is 


The  physical  position  x  =  (x^,X2>  of  the  point 


x  =  J  ft 


(A-26) 


where 


is  called  the  grid  spacing  matrix  (see  Fig.  A-l).  It  then 
follows  that  the  normalized  frequencies  p^  and  P2  can  be  identi¬ 
fied  with  the  physical  frequencies  s^  =  p^/t^  and  S2  =  P2/X2- 
Letting  s  =  (s^,S2)^  =  J  Eqs.  A-21  and  A-25  become, 

respectively, 


G(Js)  =  ^  G(ft)e"i27l<J-’^> 


(A-27 ) 


s2  S1 


G(ft ) 


=  Tlx2  f  f  G(Js)ei2n<J-'->  dSlds2  (A-28) 


s2  ”S1 


where  the  factor  1^2  in  Eq.  A-28  arises  from  the  change  in 
integration  variables.  The  quantities  5^  =  l/2tj  and  ^  =  1/2x2 
are  called  the  Nyquist  frequencies  of  the  grid. 

For  purposes  of  compatibility  with  the  definition  of 
the  Fourier  transform  in  the  continuous  case  (Eqs.  A-l  and  A-4), 
it  is  convenient  to  define  the  unnormalized  finite  Fourier 
transform  (or,  simply,  finite  Fourier  transform)  as' 


*With  this  definition  the  physical  units  of  g  and  G  coincide 
when  G  is  a  sampled  version  of  g. 


r 


f 

8 


G(s)  =  1^2  G<J§.) 


(A-29) 


Thus,  from  Eqs,  A-27  and  A-28,  the  finite  Fourier  transform  is 
given  by 


G(s)  =  A ( J )  G(Q)  e"i2n<‘I-,-:> 


and  the  inverse  finite  Fourier  transform  is 


s  ^  s  ^  a 

G(Q)  =  J  f  G(s )  ei2n< J-»->  ds1ds2 
J-s2  -s1 


(A-30) 


(A-31) 


where  in  Eq.  A-30,  A(J)  =  t^t2  is  the  determinant  of  the  spac¬ 
ing  matrix. 


A.  2  CONVOLUTION 

r 

Let  g  and  f  be  two  functions  defined  for  all  x=(x^,x2) 
The  convolution  of  g  and  f  is  a  third  function,  h,  given  by 


OB 

(x)  =  JJ  g(x' 


)  f(x-x')  dx jdx2 


(A-32) 


In  the  discrete  case,  the  convolution  of  G  and  F  is  defined  as 


H(Q)  =  G<G'  >  F(0-Q' ) 


(A-33) 


A-10 


By  changing  integration  variables  in  Eq.  A- 32  and  summation 
indices  in  Eq.  A-33  it  can  be  verified  that  the  roles  of  g  and 
f  and  those  of  G  and  F  are  interchangeable  in  the  definitions. 


A  convolution  integral  becomes  a  simple  algebraic 
product  under  the  Fourier  transformation.  To  see  this,  con¬ 
sider  Eq.  A-32.  Taking  Fourier  transforms  on  both  sides  of 
the  equality,  Eq.  A-32  becomes 

00  00 

g(x')  f (x-x' )  dxjdx£dx1dx2 

-00  -00 

(A- 34) 


Changing  the  variables  of  integration  x^  and  x2  through  the 
definition  x"=x-x' ,  Eq.  A-34  transforms  into 


00  00 

h(s)  =  ff  JJ  g(x')  f(x")  e-i2n<x’+x",s>  dxjdx^dx^dx^ 
-00  -00 


(A-35 ) 


which  can  be  written  as 


h(s)  = 


=  |yy*g(x')e’i27t<-,,->  dx|dx|j  f(x")e’i27t<-,,->dx»dx" 

(A-36) 


Thus , 


h(s)  =  g(s)  f(s) 


(A-37 ) 


Similarly,  in  the  discrete  case  the  normalized  finite 
Fourier  transform  of  H  is  given  by 


H<£)  =  G<£)  F(£) 


(A- 38) 


For  the  (unnorraalized)  finite  Fourier  transform  the  relation 
is 


H(s) 


_  1 
*  aTjT 


G(s)  F(s ) 


(A-39 ) 


As  an  example  of  a  convolution  integral,  consider  the 
flat-earth  upward  continuation  formula  of  Ref.  2.  The  poten¬ 
tial  at  height  z  over  the  point  x,  T  (x)  is  related  to  the 
surface  potential  Tq(x)  through  the  formula 


T2(x) 


z 


os 

;/ 


[||x'-x||2+z2]3/2 


To<*’ 


dx^dx^ 


(A-40) 


where  |  |x|  |  =  <x,x>  is  the  square -magnitude  of  the  vector  x. 
This  integral  can  be  recognized  as  the  convolution  of  Tq  and 
the  function  U(x)  defined  by 


U(x) 


z/2  n 

(l|x||0)V2 


(A-41 ) 


Therefore,  from  Eq.  A- 37 

T  (s)  =  U(s)  T(s)  (A-42 ) 

M  U 

At. 

It  is  instructive  to  evaluate  U(s).  More  generally, 
suppose  that  it  is  desired  to  evaluate  the  Fourier  transform 
of  a  function  r(||x||)  which  only  depends  on  the  magnitude  of 
the  position  vector  x.  Thr  Fourier  transform  r(s)  is  given  by 


r(s)  = 


r( I  I x  I  | ) 


-i2n<x , s> 
e  —  — 


dxi dx0 


(A-43 ) 


or  expanding  the  scalar  product  in  the  exponent 


r(s)  = 


00 

// 


-i2n(s,x,+s0x0) 

r( I |x| | )  e  11  1  z  dxxdx2 


(A-44 ) 


The  dependence  of  r  on  the  magnitude  of  x  suggests  the  use  of 
polar  coordinates.  Setting 


X^  -  p  cos  0 
x2  =  p  sin  0 


(A-45 ) 


Eq.  A-44  becomes 
00 

r(s)  =  f  pr(p ) 
o 


2np(s^cos0+s2sin0 ) 


d0 


dp 


(A-46) 


but  (Ref.  18) 


i27ip(s,cos0+s9sin0 ) 

e  x  c  d0  =  2n  Jo(2nsp)  (A-47) 

-n 

where  J  is  the  Bessel  function  of  the  first  kind  of  order 
o 

zero  and 


s  =  I  Is | | 


( A-48 ) 


Thus , 


r(s) 


(2nsp ) 


r(p )  dp 


( A-49 ) 


! 


I 

from  which  it  follows  that  r(s)  is  only  a  function  of  the  fre¬ 
quency  magnitude  s. 

|  The  integral  in  Eq.  A-48  is  known  in  the  mathematical 

literature  as  the  zero-order  Hankel  transform  of  the  function 
r.  From  a  table  of  Hankel  transforms  (Ref.  19),  it  follows 
that  the  Fourier  transform  of  U  as  given  by  Eq.  A-41  is 

U(s)  =  e‘2nsz  (A-50) 

Therefore,  Eq.  A-42  becomes 

Tz(s)  =  e'2;ISZ  Tq(s)  (A-51 ) 

which  is  identical  to  Eq.  A-19. 


A. 3  SPECTRAL  AND  CROSS -SPECTRAL  DENSITIES  -  DEFINITIONS 

Let  y^  and  y2  be  two  zero-mean  stochastic  processes 
defined  for  all  points  in  the  plane.  The  process  is 

said  to  be  stationary  if  its  covariance  is  only  a  function  of 
the  coordinate  differences;  i.e., 

Rv  v  U'-x")  =  E  [y«  (x  ’ )  y«(x")]  (A-52) 

y£.y£  *  * 


where  E  denotes  the  mathematical  expectation  operator.  Since 
the  roles  of  x'  and  x"  are  interchangeable  on  the  right-hand 
side  of  the  above  equation, 


*The  converse  can  also  be  shown  to  be  true:  If  r(s)  is  only  a 
function  of  ||s||,  the  inverse  Fourier  transform,  r(x)  is  only 
a  function  of  Tlx| | . 


A-14 


The  processes  and  y2  are  said  to  be  jointly  sta¬ 
tionary  if  each  of  them  is  stationary  and  their  crosscovari¬ 
ance  is  only  a  function  of  the  coordinate  differences: 


iyl ,y2^~' >  =  E*yl(-'^  y2^5")l 


(A-54) 


From  Eq.  A-54  it  follows  that 


Ry  y  <S>  =  RVl  y  <-*) 

y2’yl  yl,y2 


(A-55) 


The  power  spectral  density  (PSD)  of  y0  ,  4>  ,  is 

y£,y£ 


defined  as  the  Fourier  transform  of  its  autocovariance  function 


00 


4> 


V5>  dVx2 


(A-56) 


It  is  shown  in  Ref.  15  that 


„  (s)  >  0 


y£ ,y£  ~  ~ 


(A-57 ) 


The  cross-spectral  density  of  y,  and  y9 ,  4>  ,  is 

■L  ^  Yi  >y o 

similarly  defined  as 


yl 


,y2<->  -  ff  \.y2W  e-12"<X-'5>  d*ld*2 


(A-58) 


From  Eq.  A-55, 


Analogous  definitions  apply  to  discrete  processes. 
Let  Y^  and  Y2  be  two  discrete  zero-mean  processes  defined  for 
the  integer  coordinates  ft  =  (j,k).  The  process  Y£  is  station¬ 
ary  if 


Ry  y  te'-fl")  =  E[Yp(fl’)Y.(0")] 


(A-60) 


and  Y^  and  Y2  are  jointly  stationary  if 


1y1'Y2,<“,’-">  =  EIY1(^,)Y2(^,,)1 


(A-61) 


Their  normalized  spectral  densities  are  defined  by 

♦Y  y  <e>  =  T.  ry  y  <2>  *'12,I<-,2> 

je,Ym  Y0»Ym 


( A- 62 ) 


£ ’ ‘m 


for  £=1,2  and  m=l,2.  The  (unnormalized )  spectral  densities 


Y  (s)  =  A(J)  V] 

m 


£  ’  m 


Rv  v  (fi)e 
xl  ’  m 


-i2n< Jft , s> 


(A-63) 


A. 3.1  Spectral  Densities  and  Convolution 

Suppose  that  the  process  y2  can  be  obtained  by  con¬ 
volving  the  stationary  process  y^  with  a  known  deterministic 
function  f 


00 

y2(x')  =  f(w')  y1(x'-u»')  du>|du>^ 


(A-64 ) 


First  it  will  be  shown  that  the  covariance  of  y2 ,  E[y2(x ' )y2(x" ) 
is  only  a  function  of  the  coordinate  shift  x'-x". 


From  Eq.  A- 64, 


Ely2(x')y2(x")] 


f  (u>* ) f  (uj")E [y j  (x '  -u)'  )yj  (x"-uj" )  ] 


duj^dui^dw'j'duj” 

(A-65) 


But  since  is  stationary, 

E[y1(x'-u»')y1(x"-u)")]  =  R„  „  [  (x*  -x" )- (u»  '  -u>" )  ]  (A-66) 

i  -  -  i  yi*yi 


Thus 


00  00 


E(y2(x’)y2(x")]  =  JJ  JJ f  (u»*  ) f  (u»’*)Ry^  y ^ [  (x* -x,,)-(ui, -u>")  ]  dwjdwjdw'j'duijj 


-OB  -0D 


(A-67) 


This  equation  shows  that  the  covariance  of  y2  is  a  function  of 
x'-x".  Therefore,  y2  is  also  a  stationary  process  and  its 
covariance  is 


f(u)')f(u»") 


R  „  [x-(u)'-uj")1 
yl’yl 


dut^dw^dw'^du)^ 


(A-68) 


A  formula  for  the  spectral  density  of  y2  in  terms  of 
that  of  y^  is  derived  next.  From  Eq.  A-68, 


<t> 

y2.y2 


(s) 


If  II  If  £(->f(-"> 


Rv  v  (x-(u»'-u»")  ] 
yl’yl  -  -  ~ 


X  e-i2/t<x,s>  dw.dujt  du"dw"  dx,dx0  (A-69 ) 


Replacing  the  integration  variables  x^  and  X2  by  and  4 2 
with  4.  =  x  -  (uj  '  -w" )  and  noting  that 

e~i2n<x,s>  _  e-i27i<u»'  ,s>  e+i2n<w" , s>  e-i2n<£,s>  (a-70) 

Eq.  A- 69  becomes 

» 

♦  y  (s)  =  ff  f<"’)  e"i27l<"  ,~>  du)idw2 

-00 

00 

X  ff  f(u,»)  e  +  l2n<Ui",S> 

-00 

*  If  *y1.y1<t>  e'i2n<L,s-  dM«2 

(A-71 ) 

from  which 

4>  v  (s)  =  f<s>  **<*)  ♦v  v  <s>  <A-72) 

y2*y2  -  -  -  y1,yl 

or 

*  v  <s>  =  |f(s)|2  <t>v  v  (s)  (A-73) 

y2’y2  _  yl’yl 

which  is  the  desired  relation. 

Next,  it  is  shown  that  and  y2  are  jointly  station 
ary.  From  Eq.  A-64,  the  crosscovariance  E[y2(x' Jy^x")  ]  can 
be  written  as 

00 

E(y9(x'  )yi  (x")  ]  =  ff  f(u»)  E[y,  (x' -w)y,  (x")  Jdwndu), 


and,  because  of  the  stationarity  of  , 


00 

E(y2(x’)y1(x»)]  =  f(u>)  Ry^y^(x-u*)  du^du^  (A-75) 


Therefore,  the  crosscovariance  of  y2  and  y^  only  depends  on 
the  coordinate  difference  x.  Thus,  y^  and  y2  are  jointly  sta¬ 
tionary  and  the  crosscovariance  R  (x),  given  by  Eq.  A-75, 

y2,yl 

is  the  convolution  of  f  and  R„  „  .  It  then  also  follows  from 

yl’yl 


Eqs.  A-37  and  A-75  that  the  cross-spectral  density  4>  is 

y  2  *  y  i 

given  by 


=  f(5>  *y1,y1<5> 


(A-76) 


Similar  results  apply  to  discrete  processes.  If  the 
process  Y2  is  defined  in  terms  of  Y^  through  a  convolution 


Y2(Q.)  =  F(ft' )  Yx(fi-Q') 
Q' 


(A-77 ) 


then  the  normalized  spectral  and  cross-spectral  densities  of 
Y^  and  Y2  are  related  by 


<D 


y2,y2(^)  =  <|,y1,y1 


<£> 


(A-78) 


$ 


Y  Y  —  ^ (fi)  *^y  V 

2’*1  Yl»Yi 


<£> 


(A-79 ) 


and  the  (unnormalized)  densities  are  related  by 


*Y  Y  <®>  =  -V—  I  F<S>  | 2  $Y  Y  (s) 
2’Y2  A/(J)  Y1’y1 


(A-80 ) 


A-19 


Denoting  the  covariance  functions  in  the  primed  and  unprimed 

frames  by  R'  and  R  ,  respectively,  it  can  be  seen  that 
y » y  y » y 


-  Ry,y 


(A-84) 


The  spectral  density  of  y  expressed  in  terms  of  fre¬ 
quencies  measured  in  the  directions  of  the  primed  axes  is 


00 

*y<  ,y'  {®'  >  =  ff  Ry.,y  t '  >  dx'dx^  (A-85) 


From  Eq.  A- 84, 

00 

♦y  ’  ,y  ’  (S  ’  }  =  ff  Ry,y<0^’>  e"i2n<-'  dx^x'  (A-86) 

-» 

-1  T 

Since  6  =  6  ,  it  follows  that  <x' ,s'>  -  <0x',0s'>.  Thus, 

changing  variables  of  integration  in  Eq.  A-86  according  to  the 
transformation  A-82, 

00 

♦y .  y .  (s’ }  =  ff  Ry,y(*>  e‘i2n<-’0-’>  dxxdx2  ( A-87 ) 

-00 

or 

*’<  v» >  =  v<0^,)  (A-88 ) 

y  »y  y »y 

Equations  A-84  and  A-88  show  that,  in  general,  the 
covariance  function  and  the  PSD  of  y  depend  upon  the  specific 
directions  in  which  the  coordinate  differences  are  measured. 

A  stationary  process  whose  statistical  behavior  is  the  same  in 
all  directions  is  said  to  be  isotropic.  If  the  process  y  is 
isotropic,  the  functional  form  of  its  covariance  and  PSD  must 
be  independent  of  the  orientation  of  the  cartesian  frame  in 


which  they  are  expressed.  From  Eqs .  A-84  and  A-88  it  follows 
that  y  is  isotropic  if  for  all  possible  rotation  matrices 


Ry,y(->  =  Ry ,y(6-)  (A'89) 

and 

*  „(s)  =  (0s)  ( A-90 ) 

y  »y  y  >y 

By  changing  to  polar  coordinates,  it  can  be  shown  that  Eqs.  A-89 
and  A-90  are  possible  only  if  R  and  *  are  functions  of 

y  »y  y  »y 

| |x| |  and  s  =  | |s| |  ,  respectively  . 

An  example  of  the  covariance  of  an  isotropic  process  is 
Ry>y(x)  =  o2e  1  •  — 1  1  (1  +  p||x||  +  i  p2  I  I  x  |  |  )  ( A-91 ) 


which  is  called  an  isotropic  third-order  Markov  model  with 
2 

variance  o  and  characteristic  distance  1/p  (Ref.  14).  Using 
Eq.  A-49  and  a  table  of  Hankel  transforms,  it  is  straightforward 
to  show  that 


<D 


(S)  =  _ lOgg-V _ 

•*-  IP2  ♦  (2ns)2)7/2 


(A-92 ) 


Another  example  is  furnished  by  the  Attenuated  White  Noise 
model  (Ref.  13).  The  covariance  is  given  by 


Ry.y(s>  ■ 


8D3c6 

lllxll2  +  4  D2\'i/2 


(A-93> 


*That  these  two  conditions  are  compatible  (in  fact  one  implies 
the  other)  is  shown  by  the  derivation  of  Eq.  A-49  from  Eq.  A-43. 


where  o  is  the  variance  and  D  is  called  the  shell  depth.  The 
PSD  corresponding  to  Eq.  A-93  is 


4*  (s)  =  8tiD2ct2  e‘4*Ds  ( A-94 ) 

y  >y 

A. 3. 3  Along-track  Spectral  Densities 

The  Fourier  transform  of  a  function  of  a  single  vari¬ 
able,  g(t),  is  defined  in  a  manner  analogous  to  Eq.  A-l: 


g(f) 


-i2ntf 

e 


dt 


(A-95) 


The  inverse  one-dimensional  Fourier  transform  is 


g(t) 


00 

J  g(f)  ei2ntf  df 
-00 


(A-96 ) 


In  Eqs.  A-95  and  A-96,  f  is  interpreted  as  frequency.  The 
units  of  f  are  the  inverse  of  the  units  of  t. 


Spectral  densities  of  one-dimensional  processes  are 
denoted  by  the  letter  S  and  covariances  by  the  letter  C.  If 
z(t)  is  a  one-dimensional  process,  its  power  spectral  density 
is 


S 


z  ,z 


(f) 


/  Cz,z(t>  e'12"t£  dt 
-00 


(A-97 ) 


The  (one-dimensional)  covariance  C  _(t)  can  be  recovered  from 

the  PSD  S  „(t)  by  the  inverse  Fourier  transform  (Eq.  A-96). 
2)2 


Let  y(x)  be  a  stationary  process  on  the  plane  with 

covariance  R  (x)  and  PSD  ♦  (s).  Suppose  an  observer  is 

y  »y  y  >y 

traveling  on  the  plane  at  constant  velocity  V  on  a  straight- 
line  track  parallel  to  the  xj^  axis  of  Fig.  A-2  and  is  con¬ 
tinuously  recording  (without  errors)  the  values  of  the  process 
y  as  a  function  of  time.  Denote  by  z(t)  the  recorded  data. 

It  is  desired  to  relate  the  covariance  function  and  the  PSD  of 
z(t)  to  the  covariance  and  PSD  of  the  two-dimensional  process 
y(x) . 


First,  consider  the  covariance  C  (t).  When  the 
observer  records  an  interval  of  time  t,  the  two  points  on  the 
plane  corresponding  to  the  times  at  the  beginning  and  at  the 
end  of  the  time  interval  are  separated  by  a  vector  difference 


x  = 


cl\  /Vt  cos  ®\ 
^2  J  1  Vt  sin  0  J 


(A-98) 


In  the  primed  coordinate  system  this  difference  is 


Vt1 

0 


(A-99 ) 


Therefore  C  (t)  can  be  expressed  as 

£  )  Z 


Cz,z<C>  *  Ry,y  KVt,0)‘l 


(A-100) 


Next,  consider  the  PSD  S  (t).  Since  from  Eq.  A-85 

Z  |  2 


Oft 

R;,y<X’>  =  ff+y,  y<S’>  ei2n<-'  >  dsjds 


(A-101 ) 


A- 24 


Setting  x'  =  (Vt,0)  ,  it  follows  from  Eq.  A-100  that 


Cz,z(t> 


00 

=  /  /  $ '  (s 

JJ  y>y ~ 


,)  ei2nVtsJ  ds,ds. 


(A-102 ) 


where  sj^  is  measured  in  the  direction  of  the  x|  axis.  Substi¬ 
tuting  integration  with  respect  to  sj^  by  integration  with  re¬ 
spect  to  f=Vsj,  Eq.  A-102  becomes 


cz,z<c>  * 


00  I  00 

fy  f 


,((f/V,s£)]  ds’  (  ei2ntf  df 


(A-103 ) 


The  right-hand  side  of  Eq.  A-103  is  readily  recognized  as  the 
inverse  (one-dimensional)  Fourier  transform  of  the  quantity  in 
braces.  Therefore, 


Sz,z(t> 


00 

=  V  f  S,y(f/V'S2>  ds2 


(A-104) 


which  is  the  desired  relation. 


If  the  observer  records  the  data  as  a  function  of 
distance  rather  than  time,  the  covariance,  C  (x.'  ) ,  and  the 

y *y  i 

PSD,  S  (s,'),  of  the  resulting  record  are  called  along-track 

y »y 

covariance  and  PSD  of  the  process  y.  Formulas  for  these  quan¬ 
tities  can  be  obtained  by  setting  V  =  1,  t  =  xj  and  f  =  sj  in 
Eqs.  A-100  and  A-104: 


Cy,y(xi)  =  R;,yf<xi’0>T) 


( A- 105 ) 


A-25 


( A- 106 ) 


OB 


,y<si> 


I 


V(-' 

>y 


)  ds/ 


In  words,  the  along- track  covariance  is  the  two-dimensional 
covariance  with  the  shift  in  the  cross-track  direction  set 
equal  to  zero,  and  the  along- track  PSD  is  the  integral  of  the 
two-dimensional  PSD  over  the  cross- track  frequencies.  Note 
that  the  time-covariance  and  PSD  are  related  to  the  along- 
track  covariance  and  PSD  through  the  simple  relations 

Cz,z(t>  =  Cy.y(Vt)  u'107) 

sz,z(c>  *  V  Sy,y<£/V)  (A'108) 

In  general,  C  and  S„  depend  upon  the  direction  of 

y  *y  y  >y 

the  track..  However,  if  the  process  y  is  isotropic,  Eqs.  A-89 
and  A-90  imply  that  the  functional  forms  of  the  along- track 
covariance  and  PSD  are  independent  of  the  direction  of  the 
track. 


APPENDIX  B 

FLAT- EARTH  FREQUENCY -DOMAIN  RELATIONS 


This  appendix  presents  detailed  derivations  of  the 
flat-earth  frequency-domain  relations  introduced  in  Chapter  2. 
Relations  between  geodetic  quantities  and  between  their  sta¬ 
tistics  are  analyzed  in  Sections  B.l  and  B.2  of  this  appendix, 
respectively . 


B.l  GEODETIC  RELATIONS 


The  anomalous  potential  is  a  harmonic  function  at  all 
points  external  to  the  earth's  surface.  If  the  earth  is  taken 
to  be  a  sphere  of  radius  R,  the  anomalous  potential,  T(P) ,  at 
a  point  P  located  at  a  distance  r  >  R  from  the  center  of  the 
earth  is  given  in  terms  of  the  potential  at  the  earth's  sur¬ 
face,  Tq,  by  Poisson's  integral  formula 


T(P)  = 


_  1 


4nR 


/ 


( r2-R2 ) 


T  do 
o 


(B-l) 


where  the  integration  is  taken  over  the  surface  of  the  earth 
and  where  £  represents  the  distance  from  the  point  P  to  the 
surface  element  do. 


The  flat-earth  upward  continuation  formula  (Ref.  2) 
is  an  approximation  to  Eq.  B-l  which  corresponds  to  the  solu¬ 
tion  of  Laplace's  equation  with  a  boundary  condition  on  a  plane. 

The  potential  T  at  height  z  =  r-R  over  the  point  x  on  the 

£ 

plane  is  given  by 


B-l 


dx|dx2 


(B-2 ) 


Tz< 


aJ-.  ff  V5') 

-  (Ilx'-Xll2  ♦  .-.2I3 


T  T 

where  x  =  (x^.x^  ,  x'  =  (xj,x2>  and 

||x'-x||  =  l(xj-Xl)2  +  (x£-x2)2]1/2 


(B-3) 


The  approximation  is  valid  for  z<<R  and  | |x| |<<R. 
For  reasons  that  will  become  clear  in  the  sequel,  it  is  con¬ 
venient  to  identify  the  directions  of  the  x^  and  x2  axes  with 
the  east  and  north  directions,  respectively,  at  the  point  cor¬ 
responding  to  the  origin  of  the  plane  of  the  flat-earth 
approximation. 


The  usefulness  of  the  approximation  lies  in  the  fact 

that  the  integral  in  Eq.  B-2  is  a  two-dimensional  convolution. 

The  potential  at  height  z,  T  ,  is  the  convolution  of  the  poten 

tial  at  the  surface,  T  ,  and  the  function  U  defined  by 

o  z 


Uz(x) 


(llxll 


z/2n 


(B-4) 


The  Fourier  transform  of  a  convolution  integral  is 
equal  to  the  product  of  the  Fourier  transforms  of  the  func¬ 
tions  being  convolved  (see  Appendix  A,  Eq.  A-36).  Denoting 
Fourier  transforms  by  a  superscript  circumflex,  the  frequency- 
domain  relation  corresponding  to  Eq.  B-2  can  be  written  as 

Tz(s)  =  Uz(s)To(s)  ( B-5 ) 

T 

where  s  =  (s^,s2)  and  s^  and  s2  are  spatial  frequencies  meas¬ 
ured  in  the  east  and  north  directions,  respectively. 


The  Fourier  transform  of  U  is  given  by 

z 

o& 

Uz(s)  =  ff  U2(x>  ei2"<5.S>  dXidX2 

-00 


(B-6) 


where 

<x,s>  =  +  X2S2  (B-7) 

is  the  scalar  product  of  the  vectors  x  and  s.  This  integral 
is  evaluated  in  closed  form  in  Appendix  A.  The  result  is 

Uz(s)  =  e"2nsz  (B-8) 

where  s  =  | |s| | .  Consequently,  from  Eq.  B-5,  the  Fourier  trans¬ 
form  of  the  potential  at  height  2  is  related  to  that  at  the 
surface  by  the  simple  formula 

Tz(s)  =  e~2nsz  Tc(s)  (B-9) 

•  2n 

The  function  e  is  called  the  transfer  function 

from  potential  at  the  surface  to  potential  at  height  z.  Note 
that  this  transfer  function  is  a  simple  attenuation  factor 
whose  value  decreases  exponentially  with  frequency  and  height. 

It  is  shown  in  Appendix  A  that  Fourier  transforms  of 
partial  derivatives  of  a  function  with  respect  to  the  east  and 
north  coordinates  can  be  obtained  by  multiplying  the  Fourier 
transform  of  the  function  by  the  factors  i2rxs^  and  i2ns2»  re¬ 
spectively.  Thus,  for  example,  the  Fourier  transform  of  the 
east  and  north  components  of  the  gravity  disturbance  vector  at 
height  2,  f  and  r  ,  are  given  by  the  formulas 

X-i  Xo 


B-3 


(B-10) 


r  <s)  =  i27ts1  T2(s) 

f  (s)  =  i27is2  T2(s)  ( B- 11 ) 


Fourier  transforms  of  partial  derivatives  of  the  poten 
tial  with  respect  to  height  can  be  related  to  the  Fourier  trans 
form  of  the  potential  through  equally  simple  formulas.  For 
example,  consider  the  vertical  component  of  the  gravity  dis¬ 


turbance  vector,  T, 


Its  Fourier  transform  is  defined  by 


r 

z 


(s) 


00 

jf  rz<x)  e_i2n<*.s>  dXi<iX2 

-00 


(B-12) 


Thus , 


r 


z 


(s) 


ff  h  V*>  *1*2 

-00 


a_ 
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-00 


-i2n<x,s> 
e  —  — 


dx^dx2 


3_ 

3z 


T 


z 


(s) 


( B  - 1 3  ) 


Therefore,  from  Eq.  B-9,  it  follows  that 

rz(s)  =  -2tis  Tz(s)  (B-1A) 

Using  Eq.  B-9  it  is  possible  to  express  T  ,  r  and 

X1  x2 

V z  in  terras  of  the  Fourier  transform  of  the  anomalous  surface 
potential  as 


B-4 


<  B  - 1 5 ) 


fXi<s)  =  i2nSle"2nsz  T0(s) 
rv  (s>  =  i2ns9e"2nsz  i  (s)  (B-16) 

^2  4  O  — 

f  <s)  =  -2nse“2nsz  T(s)  (B-17) 

4  O 


The  functions  multiplying  TQ  in  Eqs.  B-15,  B-16,  and  B-17  are 
called  transfer  functions  from  anomalous  surface  potential  to 
the  east,  north,  and  vertical  components  of  the  gravity  dis¬ 
turbance  vector  at  height  z,  respectively. 

Table  B-l  presents  the  transfer  functions  from  anoma¬ 
lous  surface  potential  to  all  other  quantities  of  interest  in 
this  report.  This  table  is  identical  to  Table  2.1-1  and  is 
repeated  here  for  convenience.  To  obtain  the  Fourier  trans¬ 
form  of  any  of  the  quantities  in  the  first  column  of  the  table, 
the  Fourier  transform  of  the  anomalous  potential  at  tb',  surface 
is  multiplied  by  the  transfer  function  given  in  the  last  column 
of  the  table.  All  transfer  functions  in  Table  B-l  are  obtained 
by  replacing  the  upward  continuation  operation  and  partial 
derivatives  with  respect  to  the  east,  north,  and  vertical  co- 
ordinates  by  products  of  e  ,  i2ns^,  i2ns2  and  -2ns, 

respectively. 


Table  B-l  can  also  be  used  to  relate  any  pair  of  quan¬ 
tities  listed  in  the  first  column.  For  example,  the  frequency- 
domain  relation  equivalent  to  the  formula  of  Vening  Meinesz 
expressing  the  north  component  of  the  deflection  of  the  verti¬ 


cal  in  terms  of  the  gravity  anomaly  is  obtained  by  combining 
the  transfer  functions  in  the  third  and  fourth  rows  of  the 
table  to  yield 


fSKf'-’Z-'Si  -\ 
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TABLE  B-l 

TRANSFER  FUNCTIONS  FROM  ANOMALOUS  SURFACE  POTENTIAL 


QUANTITY 

SYMBOL 

RELATION  TO  ANOMALOUS 
POTENTIAL 

TRANSFER  FUNCTION 
FROM  ANOMALOUS  * 
SURFACE  POTENTIAL 

Anomalous  Potential  at  Height  z 

Tz 

Tz 

e*2nsz 

Undulation  of  the  Geoid 

N 

V«o 

^*0 

East  Deflection  of  the  Vertical 

n 

-OT0/9Xl)/go 

-i2nsi/g0 

North  Deflection  of  the  Vertical 

t 

-<8T0/3x2)/g0 

-i2ns2/go 

Gravity  Anomaly 

Ag 

-OV3z)lz=0  -  2Tq/R 

2ns  -  2/R 

East  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

r"i 

*V8X1 

{■}-*  -2nsz 

lzns je 

North  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

s 

3Tz/8x2 

i2ns2e-2ns2 

Vertical  Component  of  the  Gravity 
Disturbance  Vector  at  Height  z 

rz 

8Tz/8z 

-2nse~2nsz 

Gravity  Disturbance  at  Height  z 

6g 

-9Tz/8z 

2nse*2”SZ 

East-East  Gradient  at  Height  z 

r 

xixl 

82Tz/3x^ 

,_2  2  -2nsz 
-**n  s^e 

North-North  Gradient  at  Height  z 

fx2x2 

82Tz/8x2 

-An2s2e"2nsz 

Vertical-Vertical  Gradient  at 

Height  z 

rzz 

82Tz/3z2 

An2s2e*2nsz 

East-North  Gradient  at  Height  z 

Pxlx2 

82Tz/8x18x2 

/»2e  e  -2ns* 

-4n  s^S2e 

East-Vertical  Gradient  at  Height  z 

^XjZ 

82Tz/8x18z 

./_2_  - sz 

-l An  Sj se 

North-Vertical  Gradient  at  Height  z 

rXjZ 

82Tz/8x28z 

• y  2  „“2nsz 
-i bn  S2se 

*s  =  <Sj+S2>1/2 


£<  S) 


-i2ns2/gQ 
2ns -2/ft"’ 


Ag(s) 


(B-18) 


where 

9.798 


g  is  the  mean  value  of  gravity  over  the  earth  (g  - 
o  2  o 

ra/secz ) . 


B-6 


B.2 


STATISTICAL  RELATIONS 


A  zero-mean  vector  random  process  defined  on  the  plane, 
w(x)  =  Wp(x)]  ,  is  said  to  be  stationary  if  its 

covariance  matrix,  R  ,  does  not  depend  on  the  specific  values 

w  y  W 

of  the  coordinates  but”,  instead,  is  a  function  of  the  differ- 

■k 

ences  between  coordinates;  i.e., 

R  „  (x’-x")  =  E[w(x ' )  wT(x")]  (B-19) 

The  k-th  diagonal  element  of  R  is  called  the  auto- 

W  |  W 

covariance  of  the  (scalar)  random  process- corresponding  to  the 
k-th  component  of  w.  The  (k,j)  element  of  R  ,  with  k/j ,  is 

W  f  W 

called  the  crosscovariance  between  the  k-th  and  the  j-th  compo¬ 
nents  of  w.  Autocovariances  are  real,  even  and  positive  semi- 
definite  functions  (Ref.  3).  Crosscovariances  are  real  func¬ 
tions  but,  in  general,  they  are  neither  even  nor  positive 
semidef inite .  The  (k,j)  and  (j,k)  crosscovariances  are  identi¬ 
cal  except  for  a  sign  difference  in  their  arguments. 

T 

Similarly,  let  u(x)  =  (u^(x),...,  u^(x)]  be  another 

stationary  vector  random  process  on  the  plane.  The  processes 

w  and  u  are  said  to  be  jointly  stationary  if  their  covariance 

matrix,  R  ,  is  a  function  of  the  coordinate  differences; 
w ,  u 

i.e., 


Rw  „  (x’-x")  =  Elw(x')  uT(x"))  (B-20 ) 

The  Fourier  transform  of  R  is  the  spectral  density 

W  5  Vf 

matrix  of  the  process  w.  It  is  customary  to  denote  spectral 

densities  by  the  Greek  letter  <t> .  Thus,  #  =  R  and 

W  |  w  w ,  w 


*E  stands  for  the  mathematical  expectation  operator. 
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4> 

w,w 


(s) 


ff  K.*  W*'i2n<*-S-  “V*2 


(B-21) 


The  diagonal  and  off-diagonal  elements  of  4>w^w  are  called  power 
spectral  and  cross-spectral  densities,  respectively.  A  power 
spectral  density  (PSD)  is  the  average  density  of  the  random 
measure  corresponding  to  the  spectral  representation  of  a  proc¬ 
ess  (Ref.  3).  Physically,  a  PSD  can  be  interpreted  as  the 
average  distribution  of  the  energy  density  in  a  process  as  a 
function  of  frequency.  Power  spectral  densities  are  real  and 
non-negative  functions.  A  cross-spectral  density  is  a  complex¬ 
valued  function  that  indicates  the  magnitude  and  phase  of  the 
correlation  between  the  spectral  representations  of  two  proc¬ 
esses.  The  (k , j )  and  (j,k)  cross-spectral  densities  are  com¬ 
plex  conjugates  of  one  another. 


The  cross-spectral  density  matrix  of  two  jointly  sta 

tionary  processes  w  and  u,  $  ,  is  defined  in  a  similar 

zL 

fashion.  It  is  the  Fourier  transform  of  their  covariance, 
i.e.  , 


(x)e*i^n<-’2>  dx.dx. 


(B-22) 


Let  y1(x)  be  a  stationary  scalar  random  process.  It 
is  shown  in  Appendix  A  that  if  y2(x)  is  related  to  y^(x)  by  a 
transfer  function  Q(s)  in  the  form 


y2(s)  =  Q(s)  y -l ( s ) 


(B-23) 


then  the  processes  y^  and  y2  are  jointly  stationary.  Moreover, 
their  power  spectral  densities  are  related  by 


(B-24) 


These  facts  entail  a  result  of  fundamental  importance 
when  the  gravity  field  on  the  surface  of  the  earth  is  seen  as 
a  realization  of  a  stationary  random  process.  Since  all  quan¬ 
tities  in  Table  B-l  are  interrelated  by  transfer  functions, 
the  specification  of  the  covariance  function  for  any  one  of 
them  determines  completely  the  covariance  structure  for  all  of 
them. 


It  is  convenient  to  develop  a  matrix  notation  for 
expressing  the  cross-spectral  density  matrix  of  the  vector 
processes  w  and  u  in  terms  of  the  PSD  of  the  anomalous  surface 
potential  when  the  components  of  w  and  u  are  any  field-related 
quantities.  To  this  end,  suppose  that  the  k-th  component  of  w 
is  related  to  the  anomalous  surface  potential,  T  ,  through  a 
transfer  function,  G^(s) ;  i.e., 

w  =  G  Tq  (B-26) 


with 


/Gi(s)\ 

Gf?> 

\gH 


(B-27 ) 
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and  that  the  j-th  component  of  u  is  related  to  TQ  through  the 
transfer  function  Fj(s);  i.e., 


u  =  F  T 
—  —  o 


(B-28) 


with 


F1(s) 

Fj(s) 

F(s)  =  f 


rq(s> 


(B-29) 


Consider  the  (k,j)  entry,  *  ,  of  the  cross-spectral 

Wk  ’ u  j 

density  matrix  <t>  .  From  Eqs.  B-26  and  B-28,  the  transfer 

w,u 

function  from  Uj  to  w^  is  seen  to  be  the  quotient  G^/F^ .  Thus, 
from  Eq.  B-25 , 


$  =  (G,/F.)  4> 

"k-uj  k  j  uj-uj 


(B-30) 


However,  from  Eqs.  B-24  and  B-28, 


V“j  =  IFjl  "To’To 


(B-31 ) 


Therefore,  combining  the  last  two  equations, 


*wk,Uj  =  °k  Fj  ^o^o 


(B-32 ) 


It  then  follows  that  the  cross-spectral  density  matrix,  $  „ , 

is  given  in  terms  of  the  PSD  of  the  anomalous  surface  potential, 


(B-33) 


with  the  vector  transfer  functions  G  and  F  as  in  Eqs.  B-27  and 
B-29 . 


A  similar  expression  can  be  obtained  for  the  spectral 
density  matrix  .  Replacing  u  by  w  and  F  by  G  in  Eq.  B-33, 

W  f  W 

the  following  equation  is  obtained: 


♦ 

w  ,w 


(B-34) 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  C 

MULTISENSOR  SURVEY  ERROR  ANALYSIS 


This  appendix  presents  a  complete  derivation  of  the 
expressions  for  the  statistics  of  the  errors  in  the  estimates 
of  the  gravity  field  obtained  from  multisensor  data.  The  prob¬ 
lem  of  determining  the  errors  in  the  estimation  of  the  gravity 
field  from  multisensor  survey  data  can  be  formulated  in  the 
following  manner:  It  is  desired  to  characterize  the  differ¬ 
ences  (estimation  errors) 

6w(x )  =  w(x)  -  w°(x)  (C-l ) 

between  the  true  values  of  the  process  w(x),  and  the  best  esti- 
mates  w°(x).  The  components  of  w  are  any  collection  of  field- 
related  quantities. 

The  estimates,  w°(x) ,  are  obtained  from  q  data  sets 

Hk  =  {*k(x)|xeMk}  »  k=l,2,...,q  (C-2) 

corresponding  to  the  measurements  that  constitute  the  survey. 
The  set  Ek  represents  a  collection  of  scalar  measurements  of  a 
single  data  type  in  which  all  measurement  points  form  a  rectan¬ 
gular  grid  Mk<  Later  in  this  appendix  (Section  C.3),  these  q 
data  sets  will  be  classified  according  to  the  grids  on  which 
measurements  are  taken.  As  a  result,  n  classes  will  be  formed, 
where  each  class  contains  q^  sets  of  measurements  on  a  common 
grid  Mp(0=l ,2 , . . . ,n) .  Within  the  flat-earth  approximation,  the 

*In  the  mean-square  sense. 
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measurement  geometry  adapts  to  the  description  given  above  in 
most  practical  cases.  Examples  of  data  sets,  Ek,  are 

•  Gravimeter  readings  collected  at  equal 
time  intervals  by  a  survey  ship  travel¬ 
ing  at  constant  speed  along  equally 
spaced  East-West  tracks 

•  Satellite  radar-altimeter  data  on  equally 
spaced  ground  tracks.  One  set,  Z,  ,  for 

kl 

ascending  passes;  another  set,  E.  ,  for 

k2 

descending  passes. 

The  measurements,  4>k,  represent  a  scalar  linear  combi¬ 
nation  of  field-related  quantities,  u^,  corrupted  by  additive 
noise,  E^;  i.e. , 

i|»k(x)  =  uk(x)  +  Ek(x);  xeMk  (C-3) 

The  measurement  errors,  Ek,  k=l,2 . q,  are  taken  to  be  sta¬ 

tionary  gaussian  processes  independent  of  the  gravity  field. 

Formulas  for  the  spectral  density  matrix  of  the  resid¬ 
uals,  6w,  are  derived  in  this  appendix.  The  analysis  begins  in 
Section  C.l  by  considering  the  idealized  case  where  the  survey 
data  consist  of  a  collection  of  continuous  measurements  (i.e., 

Mk  =  {x|x  =  (xlfx2)  ,  -®  <  Xj^  <  * ,  -oo  <  x2  <  »}  ) .  A  general 

formula  for  the  spectral  density  matrix  of  the  errors  in  the 
estimates  obtained  with  the  optimal  linear  processor  is  derived. 

Section  C.2  the  effects  of  sampling  the  field  at  regular  in¬ 
tervals  are  examined  and  a  relation  between  the  spectral  density 
of  the  field  and  a  sampled  version  of  it  is  obtained.  This  rela¬ 
tion  is  used  in  Section  C.3  to  obtain  a  practical  solution  to  the 
general  problem  formulated  above.  The  solution  is  exact  when  all 
sets  Mk  agree  with  each  other  and  is  approximate  when  this  is  not 
the  case.  In  Section  C.4,  the  results  are  specialized  to  the 


C-2 


situation  in  which  estimates  are  sought  of  spatially  averaged 
quantities .  Section  C.5  discusses  the  evaluation  of  the  average 
errors  in  a  map  of  gravity  estimates  obtained  from  multisensor 
data. 


C.l  CONTINUOUS  MEASUREMENTS  CASE 

Consider  the  case  in  which  the  data  consist  of  a  vector 

T  T 

of  measurements  £  =  (4»^»  »|»2 » •  •  •  »#p)  at  every  point  x  =  (x^,X2>  • 

The  values  of  £  are  measurements  of  a  vector  of  field-related 

quantities,  u,  corrupted  by  stationary  additive  noise,  i.e., 

y>(x)  =  u(x)  +  £(x>  (C-4) 

The  noise  process  £  is  taken  to  be  independent  of  the  gravity 
field. 


It  is  desired  to  determine 

•  The  best  linear  estimator,  w°,  of  another 
vector  of  field-related  quantities 

T 

w  =  (w^ ,W2 , . . . ,wq)  jointly  stationary 
with  u(x) 

•  The  spectral  density  matrix  of  the  esti¬ 
mation  errors,  6w  =  w  -  w°. 

These  questions  are  dealt  with  in  Subsections  C.1.1  and  C.l. 2, 
respectively. 

C.1.1  Derivation  of  the  Optimal  Smoother 

By  a  linear  estimator  it  is  meant  that  the  estimate  of 
w  at  any  point  must  be  given  by  a  suitable  linear  combination 


C-3 


I 


of  the  data  ^(x) .  The  most  general  linear  estimator  has  the 
form 


w°(Xo>  = 


00 

// 


K(Xo,x)  £  <x)  dx^dx2 


(C-5) 


where  K  is  a  qxp  matrix  function  of  x^  and  x  to  be  determined. 

The  kernel  KCx^x)  indicates  how  to  weight  the  meas¬ 
urement  ^  at  the  point  x  in  evaluating  the  estimate  of  w  at 
the  point  3^.  Because  all  the  processes  involved  are  station- 
!  ary,  for  any  x'  the  statistical  relation  between  the  measure¬ 

ment  ij)  at  the  point  x+x'  and  the  field  w  at  the  point  x^+x'  is 
identical  to  the  relation  between  the  measurement  at  x  and  the 
field  at  x  .  Therefore,  in  evaluating  the  estimate  of  w  at 

— U  — 

3^+x'  the  measurement  at  the  point  x+x’  must  be  weighted  in 
the  same  manner  in  which  £(x)  is  weighted  in  the  computation 
of  w^x^.  Thus  for  all  x' , 

K(x<)+x' ,  x+x')  =  K  (Xq,x)  (C-6) 

This  relation  implies  that  K  only  depends  on  the  difference 
between  its  arguments.  Consequently,  Eq.  C-5  can  be  written 
as  a  convolution;  i.e., 


00 

JJ "  K(x^-x)  <J^(x)  dx^dx2 


(C-7 ) 


The  error  or  residual  in  the  estimate  of  w„  (the  n-th 


component  of  w)  at  the  point  3^  is 


n 


6wn(^o)  =  wn(*o>  •  wn(^o) 


(C-8) 
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It  is  desired  to  determine  K  so  that  the  variance  of  6wn(x  ) 
is  minimized  for  every  n=l ,2 , . . . ,q. 

Let  H,  be  the  Hilbert  space  generated*  by  the  random 

variables  for  j=l,2,...,p  (Ref.  3).  Since  w°(x  )  is  a 

j  -  n  -o 

linear  combination  of  the  data,  for  each  x  ,  w  (x  )eH,  .  The 

’  —o’  n  — o  ^ 

scalar  product  between  two  elements  of  H^,  y^  and  >  is  de- 
fined  as  the  expectation  of  their  product;  i.e., 

<yl ,y2>  =  E<yly2>  (C-9) 


and  the  norm  of  each  element  is  its  standard  deviation: 


i I y 1 1  =  (<y ,y> )1/2  =  [E(y2)]1/2 


(C-10) 


The  Hilbert  space  H^  is  a  linear  subspace  of  the  Hil¬ 
bert  space  H  generated  by  the  components  of  the  processes  w  and 
From  the  Projection  Theorem  for  Hilbert  spaces  (Ref.  20), 

it  follows  that  the  element,  w°(x  ),  of  H.  that  minimizes  the 

n  — o 

standard  deviation  of  the  difference  6wn(xQ)  =  wn(xQ)  -  w°(xq) 
is  the  projection  of  wn(xQ)  on  the  subspace  H^  and  that  the 
difference  ^w^x^)  must  be  orthogonal  to  every  element  of  H^; 
i.e.. 


Et6wn(Xo)yl  =  0 


(C-ll) 


for  every  yeH^.  Since  H^  is  generated  by  the  components  of  «£, 
Eq.  C-ll  is  equivalent  to 


*The  space  H,  contains  all  linear  combinations  of  ^(x)  and  all 

limits  in  the  mean-square  sense  of  Cauchy  sequences  of  such 
linear  combinations. 
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(C-12 ) 


E(6wn<5o)  4»j(x)]  =  0 


for  all  x  and  all  j=l,2,...,p. 


F  rom  Eq .  C  -  7 , 


P  °° 

Wn<*o>  =  S  ff  kn£<^o^,)  *2(*’)  dxldx2  (C*13) 


2  =  1  -oo 


where  knJi  is  the  component  of  K  in  row  n  and  column  2.  There 
fore,  the  expectation  on  the  left  in  Eq.  C-12  can  be  written 
as 


El^wn(^o)l^j  (— }  1  =  Elwn<5o>*j<*)l 


P  ®  _ 

7^  ffknZ  (*o'x'>  <x ’  >4»j(x)  J  dxjdx^ 
2  =  1  -» 


(C-14) 


Consequently,  from  Eq.  C-12 


p  00 

!w  ,*.<5o-£>  *  £  //knj!  <5o‘5,>  R*t.K,(5''->  dxidx2 

n  J  2  =  1  -oT  *  J 


(C-15) 


Setting  r  =  x  -x  and  changing  integration  variables  through 
o  — 

the  transformation  r'  =  x'-x,  Eq.  C-15  becomes 


(r) 


dridr2 


(C-16) 
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Since  Eq.  C-12  is  valid  for  all  x,  Eq.  C-16  is  valid 
for  all  choices  of  r.  In  matrix  form,  Eq.  C-16  becomes  a  ma¬ 
trix  convolution  equation 


os 

ff  K(r-r')  Rt>1,(r')  dr^r* 
-00 


(C-17) 


Taking  Fourier  transforms  on  both  sides  of  Eq.  C-17  results  in 


=  K<*>  *4,4(5)  <C-18> 

Therefore,  the  Fourier  transform  of  the  weighting  matrix  K  in 
Eq.  C-7  is  found  to  be  given  in  terms  of  the  spectral  densities 
of  the  measurement  and  estimated  processes  by 


K 


<t> 


-1 


(C-19) 


Condition  C-ll  is  equivalent  to  the  statement  that 
the  estimation  errors  are  uncorrelated  with  the  measurements. 
Any  two  gaussian  random  variables  which  are  uncorrelated  are 
also  independent  (Ref.  21).  Thus,  if  w  and  £  are  gaussian, 
the  errors  in  the  best  linear  estimates  are  not  only  uncorre¬ 
lated  with  the  measurements  but  independent  of  them  as  well. 


In  the  preceding  derivation,  the  attention  has  been 
restricted  to  the  class  of  linear  estimators.  Next,  it  is 
shown  that  no  nonlinear  estimator  can  further  reduce  the  vari¬ 
ance  of  the  estimation  errors  if  the  process  to  be  estimated 
and  the  measurements  are  gaussian. 

Let  w'(x  )  be  any  nonlinear  estimator  of  w(x  )  ji.e., 
w *  < 3^ )  is  some  nonlinear  functional  of  the  data  £(x)]  and  let 
6w^  be  the  error  associated  with  the  use  of  w'  in  the  estima¬ 
tion  of  the  n-th  component  of  w: 


C-7 


3 


6wA(So>  =  •  wA(*o> 


(C-20) 


The  variance  of  6w^(xQ)  , 


°8W,«W  =  E  ![un<xo)-wA<xo,l2i 
n  n 


(C-21) 


can  be  written  as 


’Sw '  ,5w' 
n’  n 


*  E{f<wn<^o)'wn(So)>  +  (wn<2o)-wA(^o)>l2} 


(C-22 ) 


where  w°  is  the  best  linear  estimator  of  wn.  From  Eq.  C-8, 

=  E{(a,W  +  <wn<V  '  wA(^o>>|2}  <C-23> 


The  expression  on  the  tight-hand  side  of  Eq.  C-23  is  the  vari¬ 
ance  of  the  sum  of  two  random  variables.  The  first  random 
variable  is  6wn(xQ),  the  error  in  the  linear  estimates.  The 
second  random  variable  is  the  difference  between  the  linear 

and  nonlinear  estimates  of  w(x  ).  Since  both  estimators  are 

n  — o 

functions  of  the  data  ^  and  since  ^  and  6wn  are  independent 
because  of  the  gaussian  assumption,  it  follows  that 


2  _  2 

f6w',6w'  ~  °6w  .6w 
n  n  n  n 


♦  e{[w°(x0)  -  w-  <xo>]2}  (C-24) 


Therefore , 


_2  .  2 
ct6w' ,6w*  -  a6w  ,6w 

n  n  nr 


(C-25 ) 


n  ’  n 


proving  the  assertion. 
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C.1.2  Spectral  Density  Matrix  of  the  Residuals 

The  covariance  matrix  of  the  residuals  is  defined  as 

Elfiw(xo)6wT(xo-x) ]  =  e|Iw(xo)-w°(xo>]  [w(xo-x)-w°(xo-x)]T} 

(C- 26) 

which  can  be  rewritten  as 

E(6w(xo)6wT(xo-x) ]  =  e|[w(xo)-w°(xo) ]  wT(xq-x)  | 

-  e|[w(X0)-W°(X0)  ]  [  <W°(3C0-3C  )  ]T  ^ 

(C- 27) 

o  X 

The  second  term  on  the  right  is  E{6w(xq)  [w  (j^-x)]  }.  Since 
w°  is  a  function  of  the  measurements  «|>,  Eq.  C-ll  implies  that 
this  expectation  is  zero.  Hence,  Eq.  C-27  reduces  to 

E[6w(xo)6wT(xo-x) ]  =  E[w(xo)wT(xo-x) ]  -  E [w°(xo)wT(xo-x) ] 

(C-28) 


But 

EIw(x0)wT(x0-x>]  =  Rw>w(x)  (C-29) 

and,  from  Eq.  C-7 


E(w°(xo)w1(xo-x) ] 


00 

=  jj  K(xQ ,x ' )  E[£(x' )wt(xo-x) ]  dxjdx£ 


00 

=  ff  K(x0-x')R;fe ^Ix-^-x')]  dx[dx£ 


$ 


which  shows  that  the  residuals  are  stationary  and  that  their 
covariance  Rfiw  6w(x)  is  given  by  the  expression  in  the  right- 
hand  side  of  Eq.~C-32. 

The  spectral  density  matrix  of  the  residuals  is  the 
Fourier  transform  of  their  covariance.  From  Eq.  C-32, 


♦6w,6w(S>  =  "'w,w<S)  -  K<2>  Vw($> 


Since , 


V«<s>  = 


it  follows  from  Eq.  C-19  that 


^6w,6w  -  ^w , w  "  ^w,jj<  ^w,^ 


(C-33 ) 


(C-3A) 


(C-35) 


This  expression  for  the  spectral  density  of  the  resid¬ 
uals  is  considerably  simplified  when  the  independence  of  the 
measurements  noise  from  the  gravity  field  is  incorporated  into 
the  formulation,  and  when  the  relations  among  the  quantities 
being  measured  and  estimated  are  taken  into  account.  From  the 
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independence  of  the  field  and  measurement  errors  and  from 
Eq.  C-4  it  follows  that 


d>  =  d> 

w,j|.  w,u 


(C-36) 


and 

<t>  s  4>  +  (ji. 

u,u 

Therefore,  Eq.  C-35  can  be  written  as 


(C-37 ) 


* 


6w  ,6w 


♦ 

w,w 


♦  [♦ 
w,ul  u,u 


-1 


(C-38) 


Now,  the  estimated  variables,  w,  and  the  measured  quantities, 
u,  are  related  to  the  anomalous  surface  potential,  T  ,  through 
vector  transfer  functions  G  and  F  as  in  Eqs.  B-26  and  B-28. 
Consequently,  from  Eqs.  B-33  and  B-34, 

♦w.w  =  5  S*  \,7a  <c'39) 


F  F*  t 
o  ’  o 


(C-40 ) 


and 


*TqiTo 


(C-41  > 


Combining  Eqs.  C-38  through  C-41  the  following  expression  for 
the  spectral  density  matrix  of  the  residuals  is  obtained: 


♦ 


5w,6w  -  - 


=  G  G  { [1  -  F  (F  F 


-1 


o’V 


f]<dt  t  } 

o  ’  o 


(C-42 ) 


There  are  two  important  observations  with  respect  to 

the  above  formula.  First,  for  any  multisensor  survey  the  vector 

transfer  functions  F  and  G  are  readily  obtained  from  the  entries 

in  Table  2.1-1.  Thus,  to  evaluate  the  spectral  density  of  the 

residuals,  one  need  only  specify  the  spectral  density  matrix 

of  the  measurement  errors,  4><.  «. ,  and  the  (a  priori)  PSD  of  the 

surface  anomalous  potential,  <*T  T  .  Second,  the  spectral 

o  *  o 

density  of  the  residuals  given  by  Eq.  C-42  has  the  same  form 
as  the  expression  for  the  (a  priori)  spectral  density  of  the 
estimated  quantities,  $  ,  as  given  by  Eq.  C- 39;  i.e.,  G  G 

W  p  w 

multiplied  by  a  scalar.”  ”Thus ,  the  residuals  bear  the  same 
relations  among  themselves  as  the  original  quantities  do. 
Moreover,  if  w  =  T  ,  then  G  =  1  and,  therefore,  the  PSD  of  the 
residual  anomalous  surface  potential,  ,  is  given  by 

the  expression  in  braces  in  Eq.  C-42;  i.e.. 


*6T  ,6T  =  [1 

o  o 


-  +  ,T  ^  1-l4>T0»T0J 


and  Eq.  C-42  can  be  written  as 


(C-43) 


^fiw.fiw  =  -  -  *6T  ,6T 

—  —  o  o 


(C-44) 


An  equivalent  expression  for  the  matrix  inverse  in 
Eq.  C-43  can  be  obtained  from  the  Matrix  Inversion  Lemma 
(Ref.  22) 


Replacing  the  matrix  inverse  in  Eq.  C- 43  by  the  right-hand 
side  of  Eq.  C-34  and  simplifying  the  resulting  expression  yields 
the  following  formula  for  the  PSD  of  the  residual  surface  anoma¬ 
lous  potential: 


6T  ,6T 
o  o 


£  *£.£.  1  * 


(C-46) 


C . 2  SAMPLING  EFFECTS  -  ALIASING 

In  order  to  analyze  the  general  problem  formulated 
in  the  introduction  to  this  appendix,  it  is  necessary  to  account 
for  the  fact  that  the  survey  data  consist  of  a  discrete  collec¬ 
tion  of  measurements.  The  purpose  of  this  section  is  to  obtain 
expressions  for  the  spectral  and  cross-spectral  density  matrices 
of  sampled  (aliased)  versions  of  vector  processes  in  terms  of 
those  corresponding  to  the  underlying  continuous  processes. 

The  problem  is  more  precisely  formulated  below. 

Denote  by  U  the  discrete  vector  random  process  cor¬ 
responding  to  having  sampled  the  process  u  on  a  rectangular 

grid  M  .  The  axes  of  the  grid  M  are  parallel  to  the  axes  of 
o  o 

a  primed  reference  frame  obtained  by  rotating  the  original 
(East-North)  frame  by  an  angle  0  as  shown  in  Fig.  C-l.  The 
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R -62449a 


*2 

(NORTH) 


the  Original  Frame  (x^,X2> 

spacings  of  the  grid  Mq  are  and  in  the  xj^  and  x£  direc¬ 
tions,  respectively.  The  grid  points  are  represented  by  pairs 
of  integers  (j  ,  k  ).  The  origin,  of  the  grid  (j  =0,  k  =0)  is 

O  O  m  O  O 

located  at  the  point  r'  =  (r|,r£;  of  the  primed  frame. 

Similarly,  let  W  represent  the  process  w  sampled  on  a 

grid  M  parallel  to  and  having  the  same  grid  spacings  as  the 

C  T 

grid  Mq  but  displaced  from  it  by  a  vector  c'  =  < c  £ ,  g  ^  >  as 

measured  in  the  primed  reference  frame  (see  Fig.  C-2). 

Two  formulas  are  sought.  First ,  an  expression  for 
the  cross-spectral  density  matrix  of  the  processes  W  and  U  in 
terms  of  the  cross-spectral  density  matrix  of  w  and  u.  Second , 
an  expression  for  the  spectral  density  matrix  of  U  in  terms  of 
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Relations  Among  the  Grid  Mc ,  (jc»kc) 
the  Grid  Mq,  <  j Q , kQ ) ,  and  the  Primed 
Coordinates  (xj.x^) 


the  spectral  density  of  u.  The  second  formula  follows  from 
the  first  by  taking  the  process  w  equal  to  the  process  u  and 
by  setting  the  displacement  vector  between  the  two  grids  equal 
to  zero.  Thus,  the  analysis  below  is  restricted  to  the  first 
formula . 


Some  definitions  are  in  order.  Let  w' ,  u*  and  R* 

-  ’  -  w,u 

denote  the  values  of  the  processes  w  and  u  and  their  covariance 
as  functions  of  the  primed  coordinates.  Similarly,  let 

w  i  u 

denote  the  cross-spectral  density  matrix  of  w  and  u  as  a  func¬ 
tion  of  frequencies  measured  in  the  x|  and  directions. 

In  the  primed  reference  frame  the  physical  position, 

T 

x^,  of  a  point  with  coordinates  fi  =  (j,k)  in  the  grid  Mq  is 


x*  =  r'  +  Jfi 
— o  —  — 

where  J  is  the  spacing  matrix  given  by 


(047 ) 


(048) 


while  the  same  integer  coordinates  in  the  grid  Mc  represent 
the  point 


x '  =  c '  +  r'  +  Jfi 

c  —  —  — 

Defining  the  rotation  matrix 


(cos  6  -sin  6\ 

sin  0  cos  8 } 


(049 ) 


(050 ) 


the  corresponding  points  in  the  original  (east-north)  frame 
are 


©jfi  (051 ) 

r  +  ejfi  (052 ) 

with  r  =  0r'  and  c  =  0c'  representing  the  position  of  the  ori¬ 
gin  of  the  grid  Mq  and  the  shift  between  the  two  grids  as  meas 
ured  in  the  original  frame. 


=  r  + 
— o  — 


xc  =  £  + 


by 


The  processes  U  and  W  are  given  in  terms  of  u'  and  w' 


U(fi)  =  u '  ( r '  +  Jfi) 


(053 ) 


016 


W(fl)  =  w'(c'  +  r'  +  JO) 


(C-54 ) 


Their  covariance, 


RW,U(-)  =  UT(fi">] 


(C-55 ) 


is  easily  related  to  the  covariance  of  the  continuous  processes 
w  and  u  through  the  use  of  Eqs .  C-53  and  0  54.  The  result  is 


Ru,u(-)  =  rw,u<£’*j2> 


(C-56) 


Spectral  and  cross-spectral  densities  of  discrete 
processes  are  denoted  by  a  superimposed  tilde.  The  cross- 
spectral  density  matrix  of  W  and  U  is  defined  as  the  finite 
Fourier  transform  of  their  covariance  (see  Appendix  A, 

Eq.  A-62 )* 

♦^(s1)  =  A(J)  ^  Rw,u<Q>  e’i2n< J— >  (C-57 ) 

Q 

T 

where  s*  =  (s^.s^)  with  s|  and  s^  frequencies  measured  in  the 
xj  and  directions,  respectively.  The  relation  inverse  to 
Eq.  C-57  is 


s2  S1 

«*.,,<«>  =  /  /  iW.U<I')  e12"02’5''  (C-: 

~s2  ’si 


00  00 


stands  for  ^  and  4(J)  is  the  determinant  of  J. 

j  =  -»  k=  -oo 
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where 


s{  =  l/2x{ 
=  l/2i£ 


(C-59) 


are  the  Nyquist  frequencies  in  the  xj  and  x£  directions,  re¬ 
spectively.  The  cross-spectral  density  matrix  <J>^  ^  is  periodic 
in  sj  and  s^  with  periods  2sj^  and  2s^,  respectively. 

The  relation  between  the  cross-spectral  density  ma¬ 
trices  of  the  sampled  and  continuous  processes  is  derived  next. 
The  basic  idea  of  the  derivation  given  below  is  to  use  the 
equality  between  the  crosscovariances  of  the  continuous  and 
discrete  processes  given  by  Eq.  C-56  to  obtain  a  relation  be¬ 
tween  the  spectral  densities. 


For  the  discrete  processes  W  and  U,  the  covariance 
matrix  Ry  ^  is  given  in  terms  of  the  spectral  density  $y  ^  in 
Eq.  C-58.  For  the  continuous  processes  w  and  u,  the  analogous 
relation  is 


T 


w,u  - 


<x'  ) 


00 

-If 


♦ 


w,u  - 


(s')  e 


i2n<x  '  ,  s ' > 


dsjds^ 


(C-60) 


It  then  follows  by  combining  Eqs.  C-56  and  C-60  that 


//  *w,u<5'> 
-00 


i2n<c '+J0  ,s  ’> 
e  -  —  — 


ds|ds£ 


(C-61 ) 


Thus,  two  different  expressions  for  Ry  jj(0)  have  been 
obtained.  In  Eq.  C-58  the  integration  is  over  a  finite  range, 
while  in  Eq .  C-61  the  integration  takes  place  over  the  whole 


C-18 


aflg 


plane.  To  compare  integrands  it  is  necessary  to  reduce  both 
integrations  to  the  same  domain.  The  integral  in  Eq.  C-61  can 
be  written  as  an  infinite  sum  of  integrals  over  equal-size 
rectangles  as 

(2m+l)s£  (2A+l)s| 

RW,U<«  =  E  /  /  *w>u<5'> 

A  (2m-l)s£  (2£-l)s| 

X  ei2«<c’+JQ,s'>  ds,ds,  (C-62) 

T 

where  A  =  (£,m)  .  Changing  integration  variables  in  such  a 
way  that  all  integrals  are  taken  over  the  domain 
[  - , s i ]  x  [-S^.s^K  Eq.  C-62  reduces  to 

ij  ii 

Rw,u<a>  =  [  I 

~^2  — 

X  ei2/i<c'+JQ,s'+J  1A>  ds,ds^  (C-63) 

The  complex  exponential  can  be  rewritten  as 

i2rt<c  ' +JQ  ,s ' +J ~^A>  _  i2n<c ' ,s '+j”^A>  i2n<Jfl,s'> 

0  ““*6  -—X0 

(C-64) 

because  the  scalar  product  <Jfl,J"^A>  =  <fl,A>  is  an  integer. 
Thus,  Eq.  C-63  becomes 


=  J  J  l£K,u<Z'+J  ±)  e1 


Zn<c s ' +J  A> 


”s2  "si  ^ 


X  ei2"<J0.s’>  ds,ds. 


(C-65 ) 


From  the  uniqueness  of  the  Fourier  transformation,  it 
follows  by  comparing  Eqs.  C-58  and  C-65  that  the  desired  rela¬ 
tion  between  the  cross-spectral  densities  of  the  discrete  and 
continuous  processes  is 


Vu<5'>  =  E  ei2n<c’  ,s’+J  V 


(C-66) 


where  A  =  (£,m)1  as  before. 

The  corresponding  relations  for  the  spectral  density 
matrices  of  processes  W  and  U  are: 


?W,V<5'>  =  E  ♦i.s/s'+J'1*; 

A 

?u,u<5'>  =  E 


(C-67 ) 


(C-68) 


Since  the  processes  w  and  u  are  interconnected  by 
transfer  functions  which  have  simple  representations  when  the 
frequencies  are  measured  in  the  east  and  north  directions,  it 
is  convenient  to  rewrite  Eqs.  C-66  through  C-68  in  the  original 
reference  frame.  These  equations  become 


Vy<2>  *  E  Vu<5+e^>  ei2„<c,s+eJ-1A> 


(C-69) 
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(C- 70) 


*W,W(-}  “  S  Vw(^+0J  1-) 


$y,u(-)  ‘  4>u,u(-+0J  1-) 


(C- 71) 


where  9  is  the  rotation  matrix  of  Eq.  C-50  and  c  =  6c'  is  the 
displacement  from  the  origin  of  the  grid  Mq  to  the  origin  of 
the  grid  Mg  as  measured  in  the  east-north  frame. 

A  graphical  interpretation  of  Eqs .  C-70  and  C-71  is 
given  in  Fig.  C-3.  The  spectral  density  matrix  of  a  sampled 
version  of  a  process,  e.g.,  TT ,  is  obtained  as  an  infinite 
superposition  of  matrix  functions  of  frequency,  <t>u  u(s+0J  A), 
each  of  which  corresponds  to  the  spectral  density  matrix  of 
the  underlying  continuous  process  translated  to  the  points 
s  =  -0J-1A  in  the  frequency  domain.  These  points  are  marked 
with  solid  dots  in  Fig.  C-3.  It  is  clear  from  Fig.  C-3  that 
the  spectrum  y  is  periodic  and  that  one  period  is  repre¬ 
sented  by  any  of  the  rectangles  in  Fig.  C-3. 

A  similar  interpretation  applies  to  the  cross-spectral 

density  matrix  of  the  processes  U  and  W  (Eq.  C-66).  In  this 

case  the  cross-spectral  density,  $  (s) ,  is  multiplied  by  the 

i2n<c  s>  —  *  — 

factor  e  — which  accounts  for  the  displacement  between 
the  grids  M  and  M  .  The  product  function  is  then  translated 
to  each  of  the  points  s  =  The  sum  of  all  these  func¬ 


tions  yields 


AVERAGE  SPECTRAL  DENSITY  OF  POST- SURVEY 
GRAVITY  RESIDUALS 


A  general  solution  to  the  problem  of  characterizing 
the  errors  in  the  estimates  of  gravity  from  multisensor  survey 


C-21 


Figure  C-3  Graphical  Interpretation  of  the  Effects 
of  Sampling  on  the  Spectrum 

data  is  given  in  this  subsection.  An  exact  solution  is  found 
for  the  case  in  which  data  from  various  sensors  are  collected 
on  a  common  rectangular  grid  Mq.  A  modification  of  the  result 
ing  expression  for  the  average  spectral  density  of  the  post¬ 
survey  gravity  residuals  yields  an  approximate  formula  for  the 
case  in  which  the  data  lie  on  different  grids. 

Consider  the  situation  in  which  the  survey  data  con¬ 
sist  of  a  measurement  vector,  ^(x^,  at  every  point  3^  on  a 


grid  Mq  as  in  Fig.  C-l.  As  before,  the  measurements  correspond 
to  a  vector  of  field-related  quantities,  u,  corrupted  by  addi¬ 
tive  noise,  E;  i.e., 

=  M<Xo>  +  E(Xo>;  xQeM0  (C-72) 

Let  the  primed  coordinate  system  be  defined  as  in  the 
last  subsection,  and  let  ^  and  U  be  the  measurements  and  the 
measured  quantities  as  functions  of  the  integer  coordinates 

-o  =  (j0’ko)T  of  the  grid  Mo- 

Since  the  process  to  be  estimated,  w,  and  the  measure¬ 
ments,  ip ,  are  gaussian,  the  best  estimate  of  w,  w°,  at  an  ar¬ 
bitrary  point  x  is  given  by  some  (unknown)  linear  combination 

d 

of  the  data ,  i.e., 

w°(xa)  =  ]£  (C-73 ) 

— o 

In  contrast  with  the  continuous  measurements  case  analyzed  in 
Section  C.l  of  this  appendix,  the  weighting  coefficients 
K(xa,fto)  depend  on  the  specific  position  of  the  point  xa  on 
the  plane.  In  fact,  the  estimation  errors 

6w(xa)  =  w(xa)  -  w°(xa)  (C-74) 

are  no  longer  stationary.  This  is  not  surprising  since,  intu¬ 
itively,  it  is  to  be  expected  that  the  rms  value  of  the  esti¬ 
mation  error  gets  smaller  as  the  point  xfl  moves  closer  to  a 
point  on  the  measurement  grid  MQ.  However,  if  another  point, 
x^,  is  chosen  so  that  it  occupies  the  same  relative  position 
with  respect  to  the  grid  points  as  x  does  (see  Fig.  C-4), 

d 

then  the  weights  K(x.,fi  )  are  identical  to  the  weights  K(x  ,Q  ) 

D  O  d  O 

except  for  a  shift  in  the  integer  coordinates.  Thus,  if  W 
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represents  the  values  of  the  field  w  as  a  function  of  the  inte 
ger  coordinates  of  a  grid  Mc  parallel  to  the  grid  Mq  as  in 
Section  C.2  (see  Fig.  C-2),  the  best  estimates  of  W  can  be 


written  as 


«°<nc>  =  Z  KCac-V  *<«„> 


(C-75 ) 


T 

where  Qc  =  (jc,kc>  is  an  arbitrary  point  of  the  grid  Mc  and 
where  K(0)  is  to  be  determined. 

A  derivation  parallel  to  that  of  Section  C.l  shows 
that  formulas  similar  to  those  obtained  in  the  continuous  meas 
urements  case  apply  to  the  discrete  case  as  well.  In  particu¬ 
lar,  the  errors  in  the  estimates 


6W(flc)  =  W(QC)  -  w°(nc) 


(C- 76) 


are  stationary  and  have  an  unnormalized  spectral  density  matrix 
given  by 


6W,6W— 


(c  * ,s ' )  =  $ 


W,W 


-  <t>  ’ 

vw,u 


[♦ 


u,u 


*E  ,E  ^ 


-1 


(♦w.u’ 


(C-77) 


where  all  the  spectral  densities  in  the  right-hand  side  are 
functions  of  s'  but,  in  addition,  ^  depends  on  the  dis¬ 
placement  c'  between  the  grids  MQ  aruT  Mc  as  indicated  by 
Eq .  C-66 . 


It  is  of  interest  to  determine  the  spectral  density 

of  the  errors  on  a  grid  M  whose  absolute  position  on  the  plane 

is  chosen  at  random.  It  suffices  to  consider  displacement 

T 

vectors  c'  =  (cj.Cj)  for  which  0  <_  cj^  <  x ^  and  0  <  c£  <  x^- 
The  spectral  density  of  the  errors  in  a  random  grid  is 


<D 


6W,6W  — 


( c  '  ,  s  '  )  dc dc 


12 


(C-78) 


where  A(J)  =  x^x£  determinant  of  the  matrix  J  of 

Eq.  C-48 . 


All  spectral  densities  appearing  on  the  right-hand 
side  of  Eq.  C-77  are  independent  of  c'  except  for  the  cross- 
spectral  density  matrix  $y  ^ .  Substituting  the  right-hand 

side  of  Eq.  C-77  for  $y(£'>s')  in  Eq.  C-78  and  using 
Eq.  C-66  results  in 


EE 


A  A 


2  1 


*sb/ dc;ac' 


0  0 


(C-79 ) 


However , 


li  12 


1  /"  /■  e12n<c',J'1(A-A')>  0  lf  -  *  - 

snr  J  Jo  e  dcidc2  j !  lf  * .  * 


1  if  A  =  A 


(080) 


Therefore  Eq.  C-79  becomes 


(C-81 ) 


Replacing  $y  y(s '  ^  the  ser^es  on  t^e  right-hand  side  of 
Eq.  C-67,  Eq.  C-81  becomes 


*)  =  ?  {♦: 


♦Au  Au(S')  =  A  {♦'  (® 

oW,oW  -  A  w,w  - 


(s'+J^A) 


-  4>'  ..(s'  +J-1A) „(s')+$;  r.(s,))"1{«'.  ..(s'  +J_ *A) ]  } 


But  jj  and  4>g  ^  are  periodic;  i.e.,  for  all  A 


(C-83 ) 

^E ,E<— ' }  =  $e,E(-'+J"1-) 

(C-84) 

Consequently,  the  spectral  density  of  the  errors 
grid  can  be  written  as 

in  a  random 

*6W,6W(-')  =  *6w,6w(-'+J  1-') 

(C-85) 

where 

4>6w,6w  =  '  w , w""  ^*w , u  l$U,U  +  *E,E]  ^w.u* 

(C-86) 

Equation  C-85  admits  an  important  interpretation. 
The  errors  in  a  random  grid  can  be  viewed  as  the  sampled  ver- 

sion  of  a  continuous  process,  6w,  whose  spectral  density  is 
given  by  Eq.  C-86.  The  spectral  density  ^  represents  the 
errors  in  a  map  of  gravity  estimates  when  no  information  is 
given  concerning  the  position  of  the  measurements  from  which 
the  map  was  obtained.  Even  if  this  information  is  provided  it 
is  shown  in  Section  C.5  that  rms  values  of  the  residuals  com¬ 
puted  on  the  basis  of  Eq.  C-86  yield  the  average  rms  of  the 
map  errors. 


The  function  ^ w  given  by  Eq.  C-86  is  called  the 

average  spectral  density  ~of  the  residuals.  The  spectral  densi¬ 
ties  in  Eq.  C-86  are  all  functions  of  frequencies,  s',  measured 
in  the  directions  of  the  primed  axes.  However,  Eq.  C-86  is 
invariant  under  rotations.  Thus,  when  all  spectral  densities 


are  expressed  as  functions  of  frequencies  in  the  east  and  north 
directions,  s,  the  relation  remains  unchanged;  i.e., 


^dw.fiw  “  ^w.w  ~  ^w,u  ^U,U  +  ^E,E^  ^w.u 


(C- 87) 


Using  Eq.  C-71,  the  average  spectral  density  of  the  residuals 
can  be  written  as 


4>,  c  =  4> 

6w,6w  w,w 


*  [♦  +  ♦.  r  +  *  I"1  4>*  (C-88) 

w,ulu,u  E,E  a, a*.  w,u 


where 


*  _  Js)  =  4>  ll(s+9J'1A) 

A/0 


(C-89 ) 


Equation  C-88  can  be  compared  with  the  analogous  ex¬ 
pression  for  the  spectral  density  of  the  residuals  in  the  con¬ 
tinuous  measurements  case  given  by  Eq.  C-38.  The  two  expres¬ 
sions  are  identical  when  the  spectral  density  of  the  errors  in 
the  continuous  measurements  case,  4>«.  «.  is  identified  with 


=  +  "a. 5 


(C-90 ) 


The  first  term  on  the  right-hand  side  of  Eq.  C-90  is 
the  finite  Fourier  transform  of  the  covariance  of  the  measure¬ 
ment  errors.  The  second  term  corresponds  to  the  sampling  (or 
aliasing)  errors.  Note  that  the  aliasing  errors  appear  in 
Eq.  C-88  as  if  they  were  indeoendent  of  the  gravity  field. 


This,  of  course,  is  not  the  case.  The  aliasing  errors  at  fre¬ 
quency  s  are  correlated  with  the  components  of  the  field  at 
frequency  s+0J  ^A  for  any  integer  vector  A  /  0.  However,  the 
aliasing  errors  at  frequency  s  are,  indeed,  independent  of  the 


components  of  the  field  at  that  same  frequency.  This  is  because 
the  spectral  representation  of  a  process  is  an  orthogonal  de¬ 
composition.  Thus,  spectral  components  at  different  frequencies 


are  independent  (Ref.  3).  Since  the  aliasing  errors  at  fre¬ 
quency  s  arise  from  field  components  at  frequencies  s+0J  ^A 
for  A  /  0,  aliasing  errors  and  field  components  at  the  same 
frequency  are  independent. 

Following  the  same  procedure  used  in  Section  C.l,  it 
can  be  shown  that  the  average  spectral  density  of  the  residuals 
is  given  by 


^fiw.fiw  -  -  *6To,6To 


(C-91 ) 


where  G  is  the  vector  transfer  function  from  anomalous  surface 


potential  to  the  quantities  being  estimated  and  where 


6T  ,6T 
o  o 


-  *£..£  1  * 


(C-92 ) 


is  the  average  power  spectral  density  of  the  residual  anoma¬ 
lous  surface  potential.  In  Eq.  C-92,  F  is  the  vector  transfer 
function  from  anomalous  surface  potential  to  the  quantities 

being  measured  and  T  is  the  a  priori  PSD  of  the  anomalous 

o  ’  o 

surface  potential. 

The  error  spectral  density  matrix  ♦  t  in  Eq.  C-92  is 
given  by  Eq.  C-90.  The  aliasing  contribution,  <l>  ,  is  ex- 

3)3 

pressed  in  terms  of  the  spectral  density  of  the  measured  quan¬ 
tities  in  Eq.  C-89.  For  computational  purposes,  it  is  conven¬ 
ient  to  obtain  it  directly  from  the  spectral  density  of  the 

unsurveyed  anomalous  surface  potential,  4>T  T  .  Using  the 

1o ’  lo 

results  of  Appendix  B.2,  Eq.  C-89  can  be  rewritten  as 

*  (s)  =  2  F(s+ej_1A)  F*(s+ej_1A)  T  (s+ej^A) 

A?0  ^ o  ’  ^ o 

(C-93 ) 

C-29 
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A  discussion  on  practical  approximations  to  this  infinite  sum 
is  given  next. 

Equation  C-93  indicates  that  the  value  of  the  spectral 

density  matrix  of  the  aliasing  errors  at  frequency  s  is  given 

by  the  infinite  superposition  of  translates  of  the  function 

—  1 

F  F  <t>T  T  to  the  points  s  =  -  0J  A  with  A/0.  This  is  illus- 
o  ’  o  ^ 

trated  in  Fig.  C-5.  The  function  F  F  <J>T  T  is  centered  at 

io*  Ao 

each  of  the  points  denoted  by  a  solid  dot  in  Fig.  C-5.  The 

aliasing  spectral  density  at  any  point  s  is  obtained  by  adding 

up  the  contribution  of  all  the  translates  of  F  F  <J>T  T  . 

1  ^  ^ 
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Figure  C-5  Graphical  Interpretation  of  Eq.  C-93 


Since  gravity  spectra  exhibit  a  rapid  decay  with  fre¬ 
quency  (Refs.  23  and  24),  in  all  practical  situations  there  is 
negligible  contribution  to  the  sum  in  Eq.  C-93  from  terms  other 

which  are  centered  at  those 


than  the  translates  of  F  F  <t>T  T 
_i  1o'io 

points  sQ  =  -©J^A  in  the  neighborhood  of  s.  In  general,  it 


is  sufficient  to  take  the  terms  corresponding  to  the  values  of 
A  belonging  to  the  set 


-1, 


Q  =  (A | A  f  0,  | | s+0J  A | |  <  2  | |s' II) 


(C-94 ) 


where  s'  =  (sj.s^)  is  the  vector  of  Nyquist  frequencies  asso¬ 
ciated  with  the  survey.  Thus,  Eq.  C-93  is  approximated  by 


(s)  =  F(s+0J  h)  F  (s+0J  ^A )  <I»T  t  (s+0J  xA) 


a ,  a 


AeQ 


To*To 


(C-95 ) 


with  Q  as  in  Eq.  C-94. 


Note  that  the  product  F  F  of  the  vector  transfer 
functions  from  surface  anomalous  potential  to  measured  quanti¬ 
ties  appearing  in  Eq.  C-95  yields  a  full  matrix.  Thus,  the 


spectral  density  matrix  of  the  aliasing  errors,  <l>  ,  is  full 

cl  f  cl 


This  means  that  aliasing  errors  in  different  measurement  types 
are  correlated.  The  reason  for  this  correlation  is  that  these 
errors  originate  from  the  same  spectral  components  of  the  field. 
This  point  is  of  importance  in  the  approximation  of  the  expres¬ 
sion  for  the  residual  spectral  density  given  below  for  the 
case  in  which  different  data  types  lie  on  nonoverlapping  grids. 
This  case  is  considered  next.  An  expression  similar  to  Eq.  C-92 
for  the  residual  anomalous  surface  potential  is  given  and  the 
approximations  involved  in  its  use  are  discussed. 
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First,  classify  the  q  data  sets  Z^  of  Eq.  02  accord¬ 
ing  to  the  grids  on  which  measurements  are  taken.  In  this 
fashion,  n  classes  result,  each  of  them  containing  qQ  sets  of 

p 

measurements  on  a  common  grid  Mp  with  0=1 ,2 , . . . ,n.  Thus, 


Q  =  +  Q2  + 


(C-96) 


If  necessary,  relabel  the  data  sets  Z^  so  that  all  sets  within 
the  same  class  are  consecutively  numbered.  The  class  identi¬ 
fied  with  the  index  3  consists  of  a  vector  of  measurements, 

,  of  dimension  qp  , 


*^q1+.  .  .+qp_1+l  ^ 

4’q1+...+qp-1+2 


(C-97 ) 


at  every  point  of  the  grid  M„ .  Let  U_  and  E„  be  the  measured 

p  ~P  ~P 

quantities  and  the  errors  in  the  measurements  of  class  p ,  and 
let  Fp  be  the  vector  transfer  function  from  anomalous  surface 
potential  to  the  quantities  being  measured  on  the  grid  M„ . 

p 

Similarly,  let  6p  and  Jp  be  the  rotation  and  spacing  matrices 
of  the  grid  Mp  defined  as  in  Eqs.  C-50  and  C-48,  respectively. 

The  approximation  to  the  average  spectral  density 
matrix  of  the  residual  anomalous  surface  potential  is  given  by 
an  expression  of  the  same  form  as  Eq.  C-92: 


6T.6T 
o  o 


F  <t> 


F  +  I/O 


To’To 


(C-98) 
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where 


(C- 99) 


with  r  the  spectral  density  of  the  errors  of  class  p 

*-p  ’*-p 
given  by 


t  <»>  =  *£  E  <£>  +  a  <£> 

^P’^P  -P~P  ~P  ~P 


(C-100) 


In  Eq.  C-100,  <!>„  _  is  the  finite  Fourier  transform  of  the 

-P  ,5P 

covariance  of  the  measurement  errors  E_  and  <J>  is  the  spec- 

P  %’~P 

tral  density  of  the  aliasing  errors  of  class  p  obtained  from  a 
formula  analogous  to  Eq.  C-95: 

«>  a  <s)  =  2  F-ts+e-J^A)  F*(s+0  J"XA)  4>T  T  (s+6  J'LA) 
V~P  “  AeQD  ~P  "  P  P  “P  ~  P  P  “  To’To  '  pP 
“  P 


(C-101) 


where  Q  is  defined  for  each  class  as  in  Eq.  C-94. 


In  writing  the  error  spectral  density  matrix,  <t> 


as  in  Eq.  C-99,  it  has  been  assumed  that  the  errors  in  the 
measurements  of  class  p,  E0 ,  are  independent  of  the  errors  in 

~p 

the  measurements  of  all  other  classes.  This  assumption  is  not 
necessary  to  obtain  an  expression  for  the  average  spectral 
density  matrix  of  the  residuals  but  it  is  satisfied  by  the 
error  models  given  in  Section  2.3,  and  it  simplifies  the  re¬ 
sulting  expressions. 


Equations  C-99  and  C-100  indicate  that  aliasing  errors 

from  different  classes  are  uncorrelated.  This,  indeed,  is  the 

case  if  any  pair  of  grids,  M0  and  MD  ,  whose  axes  are  paral- 

P1  p2 

lei  have  different  spacings  along  each  axis,  because  in  that 
case  the  aliasing  errors  in  classes  p^  and  arise  from  dif¬ 
ferent  spectral  components  of  the  field. 

The  only  approximation  involved  in  the  use  of  Eqs.  C-98 

through  C-101  (aside  from  the  flat-earth  approximation)  is  that 

the  correlation  between  aliasing  errors  of  a  given  class  and 

field  measurements  of  a  different  class  has  been  neglected. 

The  energy  at  frequency  s+0oj"^A  for  A/0  causes  aliasing  er- 

~  P  p  “  -  ~ 

rors  of  frequency  s  in  class  p.  However,  information  about 

the  spectral  component  at  frequency  s+0„J”^A  can  be  obtained 

P  P 

from  measurements  in  another  class  and  used  to  correct  for  the 
aliasing  at  frequency  s  in  class  p.  It  is  this  possible  use 
of  the  data  which  is  discounted  by  neglecting  the  correlation 
between  aliasing  errors  and  measurements  in  different  classes. 

In  most  practical  cases,  this  approximation  has  little  effect 
on  the  computed  residuals  because  the  spectral  components  af¬ 
fected  by  (correctable)  aliasing  in  class  p  are  usually  much 
better  recovered  from  measurements  in  another  class. 


Since  the  vector  transfer  function  F  in  Eq.  C-98  can 
be  expressed  as 


F 


(C-102) 


where  F^  is  the  vector  transfer  function  from  anomalous  surface 
potential  to  the  measured  quantities  in  class  0  (p=l  ,2  , . . .  ,n)  , 
Eq.  C-98  can  be  rewritten  in  a  more  convenient  form  as 


<t> 


6T  ,6T 

o  ’  o 


n 


£  *4.4  h  *  1/*Ta,T 


(S=l 


(C-103) 


Each  of  the  terms  in  the  summation  in  the  denominator  of 
Eq.  C-103  arises  from  a  single  class  of  measurements. 


The  various  measurement  classes  considered  in  this 
report  are  listed  in  Section  2.2.  Note  that  by  taking  as  many 
terms  as  appropriate  in  the  summation  of  Eq.  C-103,  any  combi¬ 
nation  of  survey  alternatives  can  be  analyzed.  In  addition, 
for  any  scalar-measurement  class  p,  the  contribution  to  the 

sum  of  the  denominator  in  Eq.  C-103  is  a  product  of  scalar 

9 

quantities  which  can  be  written  as  |F„|  /$«.  «.  .  Thus,  for 

P 

all  scalar-measurement  classes,  there  is  no  need  to  carry  out 
a  matrix  inversion  to  evaluate  the  residuals. 


C.4  RESIDUALS  OF  SPATIALLY  AVERAGED  QUANTITIES 

In  this  section,  the  formulas  given  in  Section  2.4 
for  the  average  spectral  density  of  the  residuals  of  spatially 
averaged  quantities  are  derived. 
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-  • 


Consider  a  square  block  whose  sides  are  parallel 

to  the  east  and  north  directions  and  have  length  a.  Let 

T  -a 

x  =  (x^,  X2)  be  the  center  of  the  block.  The  a-mean,  y  , 

of  a  field-related  quantity  y  at  the  point  x  is  the  average 

of  y  on  this  block,  i.e.. 


C2+a/2  x^+a/2 


ya(x)  =  K*  I  I  y(x')  dx,' dxi 

*'x2-a/2  x^-a/2  1  l 


(C-104) 


T 

where  x'  =  (xj,x£)  .  It  is  customary  to  define  a-means  only 
at  the  grid  points  of  a  square  grid  whose  axes  are  oriented 
towards  the  east  and  north  directions.  However,  it  is  con¬ 
venient  to  see  Eq.  C-104  as  defining  ya  for  all  x.  The  usual 
collection  of  a-means  is,  simply,  a  sampled  version  of  ya. 

The  quantities  ya  and  y  are  related  by  a  transfer 
function.  This  can  be  seen  by  writing  Eq.  C-104  as  a  convo¬ 
lution  ,  i.e., 


00 

ya(x)  =  JJ h(x-x' )  y(x')  dx|dx£ 


(C-105 ) 


with 


h(x")  = 


1/a4  if  | |  <  a/2  and  IxJJI  <  a/2 
0  otherwise 


(C-106) 


T 

where  x"  =  (x^.xJJ)  .  Thus,  the  transfer  function  from  y  to  y' 


00 

h(s)  =  JJ h(x")  ei2n<x",s>  dx„dx» 


(C-107 ) 
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The  evaluation  of  this  integral  is  straightforward.  The  re¬ 
sult  is 


h(s) 


sinnas^  sinpas2 
/l  TlclS^ 


(C-108) 


Next,  suppose  that  the  vector  of  quantities  to  be 
estimated  from  the  data  consists  exclusively  of  a-means.  Let 
wa  be  the  vector  of  a-means  to  be  estimated.  The  vector  trans¬ 
fer  function  from  anomalous  surface  potential  to  the  estimated 

A 

vector  is  h  G  where  G  is  the  vector  transfer  function  from 
anomalous  surface  potential  to  the  vector  w  whose  a-means  yield 
wa. 


According  to  Eq.  C-91,  the  average  spectral  density 
of  the  residuals  in  the  a-means,  <t>g^a  g^a  is 


*6wa,6wa  ~  -  -  |h|  *6T  ,6T 


(C-109) 


where  4>gT  6T  is  the  average  spectral  density  of  the  residual 
o’  o 

anomalous  surface  potential  given  by  Eq.  C-103.  Since  the 
residuals  in  the  estimates  of  w  have  an  average  spectral  den¬ 
sity  given  by 


<t’6w,6w  =  -  -  4,6To»6To 


(C-110) 


the  average  spectral  density  in  the  a-means,  w  ,  is  obtained 
from  that  of  the  residuals  in  w  by 


*6wa,6wa  =  |h|  ^Sw.fiw 


(C-lll) 


C . 5  AVERAGE  MAP  ERRORS 

The  purpose  of  this  section  is  to  prove  that,  as 
claimed  in  Section  C.3  of  this  appendix,  rms  values  of  the 
residuals  computed  on  the  basis  of  Eq.  C-86  yield  the  average 
rms  of  the  errors  in  a  map  of  gravity,  even  in  the  case  in 
which  information  concerning  the  specific  location  of  the 
measurement  grid  is  available. 

Let  the  measurement  grid  Mq  be  as  in  Section  0.2  and 
suppose  it  is  desired  to  evaluate  the  rms  of  the  errors  in  the 
estimation  of  the  vector  of  field-related  quantities  w  at  a 
point  P  as  shown  in  Fig.  C-6.  It  is  clear  that  it  is  always 
possible  to  find  a  measurement  point  with  respect  to  which  P 
is  displaced  by  less  than  Tj  units  of  distance  in  the  positive 
direction  of  the  axis  and  less  than  in  the  positive  di¬ 
rection  of  the  xj,  axis.  Let  the  displacements  in  the  xj  and 
x£  directions  be  cj  and  c£  with  0  <  cj^  <  and  0  <  c£  < 

Define  another  grid,  M£ ,  parallel  to  and  having  the  same  spac- 
ings  as  the  grid  Mq  but  displaced  from  it  by  the  amounts  c| 
and  c£  as  in  Fig.  C-7.  The  point  P  belongs  to  the  grid  Mc . 

As  shown  in  Section  C.2,  the  errors  in  the  estimation 
of  w  at  the  points  belonging  to  the  grid  Mc  are  stationary  and 
their  spectral  density  matrix  gy(£’  ' )  is  given  by  Eq.  C-77. 

Their  covariance  is  the  inverse  finite  Fourier  transform  of 

*6W,6W^-' ^ *  1,e*  » 

R6W,6W(-' 

(C-112) 


J  *SW,6W* 


e-i2n<JQ,s'>  ds,ds. 
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where  s|,  and  J  are  defined  in  Section  C.2,  Q  is  the  integer 
vector  coordinate  difference  in  the  grid  Mc  and  the  argument 
c'  of  the  covariance  emphasizes  its  dependence  on  the  shift 
between  the  estimation  and  measurement  grids. 

The  variance  of  the  estimation  errors  at  any  point  on 
2 

the  grid  Mc ,  a  (c'),  is  given  by  the  diagonal  elements  of  the 
covariance  matrix  evaluated  at  an  integer  vector  coordinate 
difference  of  zero;  i.e., 

o2(c’)  =  diag  [R6w>6w(c' ,0)]  (C-113) 

Consequently,  from  Eq.  C-112 


a2(c '  ) 


diaS[$6W,6W(^’^,)1 


ds|ds£ 


(C-m) 


The  average  variance  of  the  errors  in  the  estimates 
is  defined  as 


_1 _ 

t1t2 


dc jdc£ 


(C-115 ) 


2 

Replacing  a  (c')  by  the  expression  in  the  right-hand  side  of 
Eq.  C-114  and  interchanging  integrations  with  respect  to  dis¬ 
placements  and  frequencies  results  in 


1 

-  CM 
'  W  V 

-CM 

i  m 

r  t.  t,  -i 

f1  r 

) 

2 

a  =  diag  < 

i  i 

J  J  $6W,6W(^’’i,)dcidc2 

dsids2  } 

[-S*  -s; 

-.0  0 

) 

(C-116) 
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The  quantity  in  brackets  is  readily  recognized  as  the  spectral 
density  of  the  errors  in  a  random  grid  defined  in  Eq .  C-78. 
Thus , 


2 

cr  =  diag 


j  dsi( 


(C-117 ) 


It  was  shown  in  Section  C.3  that  <t>^y  gy  can  be  writ¬ 
ten  as  the  aliased  version  of  a  continuous  process  whose  spec¬ 
tral  density,  4>^w  gw,  is  given  by  Eq.  C-86.  From  Eqs.  C-85 
and  C-117  it  follows  that  the  average  variance  of  the  errors 
in  the  estimates  is 


2 

o  =  diag 


Si 

/  ? 


*6w,6w(£,+J‘  *>  dslds2 


(C-118) 


where  A  =  (Jl.m)1  is  a  vector  of  integers.  Interchanging  the 
operations  of  summation  and  integration,  Eq.  C-118  becomes 

®2  ®i 

a2  =  diag  £  f  I  *gw  6w(i'+J  1A)  ds^s^ 

i  m  2.  -L  -  ’  - 


"s2  _si 


(C-119) 


Consider  the  integral  I(£,m)  associated  with  the  in¬ 
dices  i  and  m  in  Eq.  C-119.  This  integral  can  be  written  as 
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5j  5[ 

I(£,m)  =  J  J  Usi+2£ii's2+2mi2)T]  dsids2 


*s2  "S1 


(C-120) 


where  use  of  Eqs .  C-48  and  C-59  has  been  made.  Changing  inte¬ 
gration  variables 

(2m+l)s2  (2£+l)s| 

I(i.m)  =  j  J  ♦'6„>Sw(2’)  dsids2 

(2m-l)i£  (2£-l)sj 


(C-121) 


Consequently,  combining  Eqs.  C-119  and  C-121, 


2 

ct  =  diag 


(2m+l)s£  (2£+l)s| 

S  2  f  f  ">6w,6w<5'>  dsids2 

<2m-l)s£  (2£-1)s| 


(C-122) 

The  integrals  in  Eq.  C-122  are  taken  over  disjoint  domains  of 
the  frequency  plane  whose  union  equals  the  whole  plane.  Thus, 


00 

o2  =  diag  f  *tw, 6„<S'>  dsids2 


(C-123) 


This  last  equation  is  the  result  that  was  sought.  It  shows  that 
the  average  variance  of  the  residuals,  defined  by  Eq.  C-115, 
can  be  computed  as  the  total  power  under  the  diagonal  elements 
of  the  average  spectral  density  matrix  of  residual  gravity. 
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